Note. Boxplots display the interquartile range (IQR, center box), and the whiskers extend 1.5*IQR from the lower and upper hinge. The white point indicates the mean and the white center line indicates the median.


Data Preparation

Data Import

# workers
# initial data cleaning was done in SPSS (syntax files are available in "")
dtWorker <- list(
  raw.pre = read_spss("data/S1_Workers/processed/cleaned/MT - Pre-Measure - 06-15-2018.sav"),
  raw.post = read_spss("data/S1_Workers/processed/cleaned/MT - Post-Measure - 06-15-2018.sav"),
  raw.morning = read_spss("data/S1_Workers/processed/cleaned/MT - Morning - 06-15-2018.sav"),
  raw.afternoon = read_spss("data/S1_Workers/processed/cleaned/MT - Afternoon - 06-15-2018.sav")
)

# students
dtStudents <- list(
  raw.pre = read.csv(file = "data/S2_Students/raw/AOTS_Pre.csv", header = T, sep = ","),
  raw.post = read.csv(file = "data/S2_Students/raw/AOTS_Post.csv", header = T, sep = ","),
  raw.daily = read.csv(file = "data/S2_Students/raw/AOTS_Daily.csv", header = T, sep = ",")
)

# young medical professionals
dtMedical <- list(
  raw.pre = c(),
  raw.post = c(),
  raw.daily = c()
)

Data Cleaning & Data Exclusions

#  important names for Morning and Afternoon
names.m <- c(
  "StartDate",
  "EndDate",
  "Finished",
  "Duration__in_seconds_",
  "RecordedDate",
  "ExternalReference",
  "Meta_Operating_System",
  "Contact_dum",
  "number",
  "time",
  "duration_1",
  "dyad.group",
  "gr_size",
  "gr_type_1",
  "gr_type_2",
  "gr_type_3",
  "gr_type_4",
  "gr_type_5",
  "gr_type_6",
  "gr_type_7",
  "gr_type_8",
  "gr_type_9",
  "gr_type_10",
  "gr_type_11",
  "gr_type_12",
  "gr_type_13",
  "gr_type_14",
  "gr_type_15",
  "gr_type_16",
  "gr_type_17_TEXT",
  "gr_context_1",
  "gr_context_2",
  "gr_context_3",
  "gr_context_4",
  "gr_context_5",
  "gr_context_6",
  "gr_context_7",
  "gr_context_8",
  "gr_context_9",
  "gr_context_10",
  "gr_context_11",
  "gr_context_12",
  "gr_context_13_TEXT",
  "gr_context_14_TEXT",
  "gr_dutchness",
  "dyad_type_1",
  "dyad_type_2",
  "dyad_type_3",
  "dyad_type_4",
  "dyad_type_5",
  "dyad_type_6",
  "dyad_type_7",
  "dyad_type_8",
  "dyad_type_9",
  "dyad_type_10",
  "dyad_type_11",
  "dyad_type_12",
  "dyad_type_13",
  "dyad_type_14",
  "dyad_type_15",
  "dyad_type_16",
  "dyad_type_17_TEXT",
  "Context_1",
  "Context_2",
  "Context_3",
  "Context_4",
  "Context_5",
  "Context_6",
  "Context_7",
  "Context_8",
  "Context_9",
  "Context_10",
  "Context_11",
  "Context_12",
  "Context_13_TEXT",
  "Context_14_TEXT",
  "keyMotive",
  "keymotive_fulfillemt_1",
  "keyMotive_Dutch_1",
  "autonomy_1",
  "competence_1",
  "relatedness_self_1",
  "relatedness_other_1",
  "qualityAccidental_1",
  "qualityVoluntary_1",
  "qualityCooperative_1",
  "qualityDutchy_1",
  "quality_overall_1",
  "quality_meaning_1",
  "quality_star_1",
  "wantInt",
  "desire_type_1",
  "desire_type_2",
  "desire_type_3",
  "desire_type_4",
  "desire_type_5",
  "desire_type_6",
  "desire_type_7",
  "desire_type_8",
  "desire_type_9",
  "desire_type_10",
  "desire_type_11",
  "desire_type_12",
  "desire_type_13",
  "desire_type_14",
  "desire_type_15",
  "desire_type_16",
  "desire_type_17_TEXT",
  "desire_context_1",
  "desire_context_2",
  "desire_context_3",
  "desire_context_4",
  "desire_context_5",
  "desire_context_6",
  "desire_context_7",
  "desire_context_8",
  "desire_context_9",
  "desire_context_10",
  "desire_context_11",
  "desire_context_12",
  "desire_context_13_TEXT",
  "desire_context_14_TEXT",
  "Reason_nodesire",
  "keyMotive_noInt",
  "keyMotive_noInt_fulf_1",
  "autonomy_NoInt_1",
  "competence_NoInt_1",
  "relatedness_1_NoInt_1",
  "thermometerDutch_1",
  "thermometerDutchInt_2",
  "ExWB_1",
  "alertness1",
  "calmness1",
  "valence1",
  "alertness2",
  "calmness2",
  "valence2",
  "inNonDutch",
  "NonDutchNum",
  "NonDutchType_1",
  "NonDutchType_2",
  "NonDutchType_3",
  "NonDutchType_4",
  "NonDutchType_5",
  "NonDutchType_6",
  "NonDutchType_7",
  "NonDutchType_8",
  "NonDutchType_9",
  "NonDutchType_10",
  "NonDutchType_11",
  "NonDutchType_12",
  "NonDutchType_13",
  "NonDutchType_14",
  "NonDutchType_15_TEXT",
  "date",
  "time.0",
  "LocationLatitude",
  "LocationLongitude"
)

names.a <- c(names.m, "keyInteraction_1", "keyInteractionTime")

# Create reduced data sets for morning and afternoon
dat.mo <- dtWorker$raw.morning[, names.m]
dat.mo$daytime <- "morning"

dat.af <- dtWorker$raw.afternoon[, names.a]
dat.af$daytime <- "afternoon"

# merge morning and afternoon measurements with indicator [+ clean up]
daily.dat <- rbind.fill(dat.mo, dat.af)
daily.dat <- daily.dat[daily.dat$ExternalReference != 55951, ]
dtWorker$daily <- daily.dat
rm(dat.mo, dat.af, names.m, names.a, daily.dat)


# names for pre-measurement
names.pre <- c(
  "Finished",
  "age",
  "Gender",
  "Living",
  "roommate_1",
  "roommate_2",
  "roommate_3",
  "nationality",
  "SecondNationality",
  "timeNL_1",
  "Reason_2",
  "Reason_5",
  "Reason_7",
  "Reason_8_TEXT",
  "DutchLang",
  "occupation_1",
  "occupation_2",
  "occupation_3",
  "occupation_4",
  "occupation_7",
  "CurrentEducation_1",
  "education_level",
  "EduLang_2",
  "RUG_faculty",
  "Study.0",
  "association",
  "DutchMeetNum",
  "DutchFriends_1",
  "assimilation",
  "separation",
  "integration",
  "marginalization",
  "VIA_heritage",
  "VIA_Dutch",
  "SSAS_surrounding",
  "SSAS_privat",
  "SSAS_public",
  "autonomy",
  "relatedness",
  "competence",
  "anxiety",
  "swl",
  "alertness",
  "calmness",
  "valence",
  "date",
  "time",
  "City",
  "ZIP",
  "id"
)

# reduced data set for pre measurement
dat.pre.red <- dtWorker$raw.pre[, names.pre]

# merge with daily data [+ clean up]
df.pre <- merge(
  x = dtWorker$daily,
  y = dat.pre.red,
  by.x = "ExternalReference",
  by.y = "id",
  all = T
)
rm(names.pre)

# adjust duplicate names to fit to indicate daily or pre measurement
names(df.pre) <- gsub("[[:punct:]]x", ".daily", names(df.pre))
names(df.pre) <- gsub("[[:punct:]]y", ".pre", names(df.pre))

# names for post measurement
names.post <- c(
  "ExternalReference",
  "assimilation",
  "separation",
  "integration",
  "marginalization",
  "VIA_heritage",
  "VIA_Dutch",
  "anxiety",
  "swl",
  "rosenberg",
  "social_support",
  "stress",
  "discrimination",
  "discrimination_month",
  "NLE_1month",
  "NLE_6month",
  "NLE_12month"
)

# reduced data set for post-measurement
dat.post.red <- dtWorker$raw.post[, names.post]

# merge post measurement with pre- and daily data
df <- merge(
  x = df.pre,
  y = dat.post.red,
  by.x = "ExternalReference",
  by.y = "ExternalReference",
  all = T
)

# adjust duplicate names to indicate pre or post
names(df) <- gsub("[[:punct:]]x", ".pre", names(df))
names(df) <- gsub("[[:punct:]]y", ".post", names(df))

# add to list
dtWorker$combined <- df

# create data frame with cleaned data
df <- dtWorker$combined %>%
  filter(
    Finished.pre == 1,
    Finished.daily == 1,
    !is.na(ExternalReference)
  )

# add running number as measurement ID within participants
df$measureID <- rowidv(df, cols = c("ExternalReference"))

df <- df %>%
  mutate(
    PID = as.numeric(factor(ExternalReference)),
    # participant ID
    TID = measureID - 1,
    # time ID with t0 = 0 for meaningfull intercept interpretations
    date = substr(StartDate, 1, 10),
    # awkward way of extracting date (best converted to )
    time = substr(StartDate, 12, 19),
    # awkward way of extracting time
    daynum = as.numeric(factor(date)),
    # all days as numeric for ordering
    daycor = ifelse(
      daytime == "morning" &
        period_to_seconds(hms(time)) < period_to_seconds(hms("12:00:00")) |
        daytime == "afternoon" &
          period_to_seconds(hms(time)) < period_to_seconds(hms("19:00:00")),
      daynum - 1,
      daynum
    ),
    # correctly identify which date the questionnaire is about
    daycor.lead = sprintf("%02d", daycor),
    daytime.lt = ifelse(daytime == "morning", "a", "b"),
    # morning / afternoon to a / b
    day_time = paste(daycor.lead, daytime.lt, sep = "_"),
    # combine day id with morning / afternoon
    session = as.numeric(factor(day_time)),
    # day and time identifier as numeric id
    SubTime = chron::times(time.0),
    time.daily = as.character(time.daily),
    PPDate = as.Date(df$date.daily),
    number = replace_na(number, 0),
    NonDutchNum = replace_na(NonDutchNum, 0)
  )

dtWorker$clean <- df

# clean up
rm(df.pre, names.post, dat.post.red, dat.pre.red, df)

# Export reduced Data
# write.csv(dtWorker$clean, "data/processed/MT_clean-merged_07-05-2018.csv", row.names = F)
# save(dtWorker$clean, file = "data/processed/MT_clean-merged_07-05-2018.RData")
# our own test IDs
ownIDs <- c(
  "beautifulLionfishXXXR5rcgVBzGu8hPvOqrK8UBJBw4owvi9nfRFSFu3lMzYhE",
  "niceDogoXXXmB8JI5SFu78SF3DVof84mGUPPNUr14p2HYFTtp31a6D1OwAzM6F-K",
  "amusedQuailXXXmhuc_fpTp8vPkMwDH1BzjaH1d1kHSO1bsPEfsnaEYk4WeVBfPi",
  "juwGAbtXX0_1kmZtSVqKh3PGaHOICqUyU4iBkrT3nDsI_uifuD1gzKcZerxaM5FL"
)

# Prepare dfs for Cleaning
df.pre <- dtStudents$raw.pre %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!is.na(ended)) %>% # remove all who did not finish
  filter(!e_mail %in% .$e_mail[duplicated(.$e_mail)]) %>% # remove all who did the pre questionnaire multiple times (b/c inconsistent ratings scales)
  filter(!session %in% ownIDs) %>% # remove our own test
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

df.post <- dtStudents$raw.post %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!is.na(session)) %>% # remove own test runs
  filter(!session %in% ownIDs) %>% # remove our own test
  filter(session %in% df.pre$session) %>% # remove anyone who wasn't in the pre
  filter(!is.na(ended)) %>% # remove all who never finished
  filter(!session %in% .$session[duplicated(.$session)]) %>% # remove all duplicate sessions
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

df.daily <- dtStudents$raw.daily %>%
  mutate_all(na_if, "") %>%
  mutate_all(na_if, "NA") %>%
  filter(!session %in% ownIDs) %>% # remove our own test
  filter(session %in% df.pre$session) %>% # remove anyone who wasn't in the pre
  filter(!is.na(ended)) %>% # remove all who never finished
  mutate(session = as.character(session)) # turn factor into character strings (probably just precaution)

# merge daily with pre
dfPreDaily <- merge(
  x = df.daily,
  y = df.pre,
  by = "session",
  suffixes = c(".daily", ".pre"),
  all = F
)

# merge daily with post
dfCombined <- merge(
  x = dfPreDaily,
  y = df.post,
  by = "session",
  suffixes = c(".pre", ".post"),
  all = F
)

# add to list
dtStudents$clean <- dfCombined

# clean up workspace
rm(df.pre, df.daily, df.post, dfPreDaily, dfCombined, ownIDs)
# TBD

Calculate needed transformations

Worker

df <- dtWorker$clean

# Time and Date Variables
# remove seconds from afternoon time
df$SubTime[df$daytime == "afternoon"] <- paste0(substring(as.character(df$time.0[df$daytime == "afternoon"]), 4, 8), ":00")
df$time.daily[df$daytime == "afternoon" &
  !is.na(df$time.daily != "<NA>")] <- paste0(substring(as.character(df$time.daily[df$daytime == "afternoon" &
  !is.na(df$time.daily != "<NA>")]), 4, 8), ":00")

# Correct morning / afternoon date where survey was collected the day after to indicate the correct date that was targeted
df$PPDate[df$SubTime < "11:50:00" &
  df$daytime == "morning"] <- df$PPDate[df$SubTime < "11:50:00" &
  df$daytime == "morning"] - 1
df$PPDate[df$SubTime < "18:50:00" &
  df$daytime == "afternoon"] <- df$PPDate[df$SubTime < "18:50:00" &
  df$daytime == "afternoon"] - 1

# Mood scales
df$calmness.daily <- rowSums(df[, c("calmness1", "calmness2")], na.rm = T)
df$alertness.daily <- rowSums(df[, c("alertness1", "alertness2")], na.rm = T)
df$valence.daily <- rowSums(df[, c("valence1", "valence2")], na.rm = T)

# Need scales
df$keyMotiveFulfilled <- rowSums(df[, c("keymotive_fulfillemt_1", "keyMotive_noInt_fulf_1")], na.rm = T)
df$autonomy.daily.all <- rowSums(df[, c("autonomy_1", "autonomy_NoInt_1")], na.rm = T)
df$competence.daily.all <- rowSums(df[, c("competence_1", "competence_NoInt_1")], na.rm = T)
# cor(df$relatedness_other_1, df$relatedness_self_1,use="complete.obs")
df$relatedness.daily.all <- rowMeans(df[, c(
  "relatedness_other_1",
  "relatedness_self_1",
  "relatedness_1_NoInt_1"
)], na.rm = T)
df$relatedness_1 <- rowMeans(df[, c("relatedness_other_1", "relatedness_self_1")], na.rm = T)

# summarize by participant (check that everything is within pp might not be the case for )
between <- df %>%
  group_by(ExternalReference) %>%
  mutate(
    CtContactNL = sum(Contact_dum),
    CtContactNonNl = sum(inNonDutch),
    CtContactNLAll = sum(number),
    CtContactNonNlAll = sum(NonDutchNum),
    AvKeyNeed = mean(keyMotiveFulfilled, na.rm = T),
    AvKeyNeedInt = mean(keymotive_fulfillemt_1, na.rm = T),
    AvKeyNeedNoInt = mean(keyMotive_noInt_fulf_1, na.rm = T),
    AvAutonomy = mean(autonomy.daily.all, na.rm = T),
    AvCompetence = mean(competence.daily.all, na.rm = T),
    AvRelatedness = mean(relatedness.daily.all, na.rm = T),
    AvThermo = mean(thermometerDutch_1, na.rm = T),
    AvWB = mean(ExWB_1, na.rm = T)
  ) %>%
  ungroup() %>%
  mutate(
    CtContactNL_c = scale(CtContactNL, scale = FALSE),
    AvKeyNeedInt_c = scale(AvKeyNeedInt, scale = FALSE),
    AvKeyNeed_c = scale(AvKeyNeed, scale = FALSE),
    CtContactNL_z = scale(CtContactNL, scale = TRUE),
    AvKeyNeedInt_z = scale(AvKeyNeedInt, scale = TRUE),
    AvKeyNeed_z = scale(AvKeyNeed, scale = TRUE)
  )

warning(
  "some variable transformations (esp. _c and _z) might be across all participants (i.e., not within PP)."
)

dtWorker$full <- between
rm(df, between)

# Between participants contact frequency
workerContactFreq <- dtWorker$full %>%
  group_by(ExternalReference) %>%
  summarise(
    n = n(),
    SumContactNL = sum(Contact_dum),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(number),
    AvAttitude = mean(thermometerDutch_1, na.rm = T)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# dataframe where interaction types are recoded
workerInteractionType <- dtWorker$full %>%
  mutate(
    OutgroupInteraction = as_factor(Contact_dum),
    NonOutgroupInteraction = as_factor(inNonDutch)
  )

# save cleaned data
# save(df.btw, file = "data/processed/df.btw.RData")
# write_sav(df.btw, "data/processed/MT_clean-merged_pre-post.sav")

# export data to Mplus
# df.mplus = remove_all_labels(select(df,
#                                     PID, session,
#                                     thermometerDutch_1, inNonDutch, Contact_dum,
#                                     keyMotiveFulfilled, autonomy.daily.all, competence.daily.all, relatedness.daily.all))
# names(df.mplus)= c("PID", "session", "att", "intin", "intout", "keymot", "aut", "comp", "rel")
# mplus = df.mplus[order(df.mplus$PID, df.mplus$session),]
# mplus.intcont = mplus[mplus$intout==1,]
# prepareMplusData(mplus.intcont, "data/processed/dynamic-subset-intonly.dat")

Student

df <- dtStudents$clean

# Add ID variables
df$PID <- as.numeric(factor(df$session)) # participant ID

# order time
df$TID <-
  factor(df$date_period, levels = unique(dtStudents$raw.daily$date_period))
df$TIDnum <- as.numeric(df$TID) # get numeric TID

# check whether time ordering worked
df <- df %>%
  arrange(PID, TID) # %>%
# View()

# Interaction as Factor
df$interaction.f <-
  factor(df$Interaction,
    levels = c("no interaction", "Dutch", "Non-Dutch")
  )
df$intNL <- ifelse(df$Interaction == "Dutch", 1, 0)
df$intNonNL <- ifelse(df$Interaction == "Non-Dutch", 1, 0)

# -------------------------------------------------------------------------------------------------------------
#                                       Combine Variables
# -------------------------------------------------------------------------------------------------------------
# Relatedness
pairs.panels.new(
  df[c("RelatednessSelf", "RelatednessOther")],
  labels = c(
    "I shared information about myself.",
    "X shared information about themselves."
  )
)

df$RelatednessInteraction <-
  rowMeans(df[c("RelatednessSelf", "RelatednessOther")], na.rm = T)
df$RelatednessInteraction[df$RelatednessInteraction == "NaN"] <-
  NA
# Relatedness Overall (JANNIS NOT SURE THESE ARE CORRECT, CHANGE ROWS?; J: Changed "NaN" in df$RelatednessInteraction to NA() should work now)
df$Relatedness <-
  rowMeans(df[, c("RelatednessInteraction", "RelatednessNoInteraction")],
    na.rm =
      T
  )
# Pro-Sociality
df$ProSo <-
  rowMeans(df[, c("ProSo1", "ProSo2", "ProSo3", "ProSo4")], na.rm = T)
# Anti-Sociality
df$AntiSo <-
  rowMeans(df[, c("AntiSo1", "AntiSo2", "AntiSo3", "AntiSo4")], na.rm = T)


# -------------------------------------------------------------------------------------------------------------
#                                 Add Variables related to interaction partner
# -------------------------------------------------------------------------------------------------------------
# create function for later lapply
createIntPartDf <- function(inp) {
  # prepare the dataframe so that we can forloop over it later
  tmp <- data.frame(
    CC = as.character(inp$CC),
    NewCC = as.character(inp$NewCC),
    NewName = as.character(inp$NewName),
    NewCloseness = inp$NewCloseness,
    NewGender = inp$NewGender,
    NewEthnicity = as.character(inp$NewEthnicity),
    NewRelationship = as.character(inp$NewRelationship)
  )

  tmp$CC2 <- recode(tmp$CC, "SOMEONE ELSE" = "NA")
  tmp$CC2 <-
    ifelse(
      tmp$CC == 1 |
        tmp$CC == "SOMEONE ELSE",
      as.character(tmp$NewName),
      as.character(tmp$CC2)
    )
  # maybe add [[:space:]]\b to remove space before word boundary or ^[[:space:]] to remove space in the beginning of a string
  tmp$CC2 <- gsub("^[[:space:]]", "", tmp$CC2)
  tmp$NewName <- gsub("^[[:space:]]", "", tmp$NewName)

  # open the variables that will be filled up in the foor-loop
  tmp$closeness <- rep(NA, nrow(tmp))
  tmp$gender <- rep(NA, nrow(tmp))
  tmp$ethnicity <- rep(NA, nrow(tmp))
  tmp$relationship <- rep(NA, nrow(tmp))

  # Run the for-loop. It finds the variables related to the name of the interaction partner. If there is a repeating interaction
  # partner (i.e. CC2) it takes the value (i.e. NewCloseness) from the first interaction (i.e. NewName)
  for (i in 1:nrow(tmp)) {
    if (is.na(tmp$CC2[i])) {
      next
    } else {
      tmp$closeness[i] <-
        na.omit(tmp$NewCloseness[as.character(tmp$CC2[i]) == as.character(tmp$NewName)])[1] # find closeness where CC2 matches NewName (na.omit + [1] to get the number)
      tmp$gender[i] <-
        na.omit(tmp$NewGender[as.character(tmp$CC2[i]) == as.character(tmp$NewName)])[1] # (na.omit + [1] to get the number and not the rest of the na.omit list)
      tmp$ethnicity[i] <-
        na.omit(as.character(tmp$NewEthnicity[as.character(tmp$CC2[i]) == as.character(tmp$NewName)]))[1] # PROBLEM IS THAT THERE ARE TOO MANY NA's: Difficult to deal with
      tmp$relationship[i] <-
        na.omit(as.character(tmp$NewRelationship[as.character(tmp$CC2[i]) == as.character(tmp$NewName)]))[1]
    }
  }

  out <- tmp
  out
}

# split df per participants and run function
PP <- split(df, df$PID)
PP <- lapply(PP, createIntPartDf)
rm(createIntPartDf)

# add variables back to df
remergePP <- do.call(rbind.data.frame, PP)
colnames(remergePP) <-
  paste(colnames(remergePP), "_Calc", sep = "")
df <- cbind(df, remergePP)
rm(remergePP, PP)

# -------------------------------------------------------------------------------------------------------------
#                                 Center Relevant Variables
# -------------------------------------------------------------------------------------------------------------

df <- df %>%
  group_by(PID) %>%
  mutate(
    KeyNeedFullfillment.cm = mean(KeyNeedFullfillment, na.rm = TRUE),
    # cluster mean (mean of PP)
    KeyNeedFullfillment.cwc = KeyNeedFullfillment - KeyNeedFullfillment.cm,
    # cluster mean centered (within PP centered)
    closeness.cm = mean(closeness_Calc, na.rm = TRUE),
    closeness.cwc = closeness_Calc - closeness.cm
  ) %>%
  ungroup()

# store
dtStudents$full <- df
rm(df)

# Between participants contact frequency
studentContactFreq <- dtStudents$full %>%
  group_by(PID) %>%
  summarise(
    n = n(),
    SumContactNL = sum(InteractionDumDutch),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(ContactNum[InteractionDumDutch == 1], na.rm = TRUE),
    AvAttitude = mean(AttitudesDutch, na.rm = TRUE),
    AvQuality = mean(quality_overall, na.rm = TRUE)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# dataframe where interaction types are recoded
studentInteractionType <- dtStudents$full %>%
  mutate(
    NonDutchContact = replace_na(NonDutchContact, 2), # make second non-Dutch countable
    NonDutchContact = NonDutchContact*-1+2 # recode (yes = 1 -> 1, no = 2 -> 0)
  ) %>%
  mutate(
    OutgroupInteraction = factor(
      InteractionDumDutch,
      levels = c(0, 1),
      labels = c("No", "Yes")
    ),
    NonOutgroupInteraction = factor(
      rowSums(select(., c(InteractionDumNonDutch, NonDutchContact))), # combine the two non-Dutch Q.,
      levels = c(0, 1),
      labels = c("No", "Yes")
    )
  )

# select a subset of IDs to display in plots
studentPltIDs <-
  studentInteractionType %>%
  group_by(PID) %>%
  summarise(n = n()) %>%
  slice_max(n, n = 20) %>% # chose the 20 with the most number of measurements
  select(PID) %>%
  as.matrix %>%
  as.vector

# select a subset of IDs to display in plots (only outgroup interactions)
studentOutPltIDs <-
  studentInteractionType %>%
  filter(OutgroupInteraction == "Yes") %>%
  group_by(PID) %>%
  summarise(n = n()) %>%
  slice_max(n, n = 20) %>% # chose the 20 with the most number of measurements
  select(PID) %>%
  as.matrix %>%
  as.vector

Medical

# TBD

Worker Sample

Data Description

Still in ‘scr/workerDescriptive.R’. Needs to be merged with this document.

Contact Hypothesis

Interaction Frequency and Outgroup Attitudes

\[\begin{equation} \tag{1} r_{ContactFrequency, OutgroupAttitudes} \neq 0 \end{equation}\]

# correlation panel
pairs.panels.new(
  workerContactFreq %>% select(SumContactNL, SumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  workerContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

TL;DR: Neither the number of interactions nor the number of measurement beeps with an interaction are significantly related with the average outgroup attitudes. This is to say that within our data participants with more outgroup interactions did not have significantly more positive outgroup attitudes. This might be due to the aggregation within the participants or the small sample size of between participant data. Nonetheless, the aggregate data does not support the notion that simply having more interactions with an outgroup results in more positive outgroup attitudes.

Outgroup Attitudes by Interaction Type

\[\begin{equation} \tag{2} \mu_{i,0utgroupInteraction} > \mu_{i,IngroupInteraction} \\ Attitude = OutgroupInteraction \times NonOutgroupInteraction \end{equation}\]

Regression with interaction term

Across participants: impact of different interaction types [probably meaningless because ignores response sets between participants]

# between participants interaction type
workerContactType <- dtWorker$full %>%
  group_by(
    Contact_dum,
    inNonDutch
  ) %>%
  summarise(
    n = n(),
    AttitudeM = mean(thermometerDutch_1, na.rm = T),
    AttitudeSD = sd(thermometerDutch_1, na.rm = T),
    AttitudeSE = AttitudeSD / sqrt(n),
    AttitudeLwr = AttitudeM - 1.96 * AttitudeSE,
    AttitudeUpr = AttitudeM + 1.96 * AttitudeSE
  ) %>%
  ungroup() %>%
  mutate(InteractionType = paste(
    ifelse(Contact_dum == 1, "Out", "NoOut"),
    ifelse(inNonDutch == 1, "In", "NoIn"),
    sep = ", "
  ))

# plot bar chart
ggplot(
  workerContactType,
  aes(
    y = AttitudeM,
    x = as_factor(Contact_dum),
    fill = as_factor(inNonDutch)
  )
) +
  geom_bar(
    stat = "identity",
    color = "black",
    position = position_dodge()
  ) +
  geom_errorbar(aes(ymin = AttitudeM, ymax = AttitudeUpr),
    width = .2,
    position = position_dodge(.9)
  ) +
  labs(
    fill = "Non-Outgroup Interaction",
    x = "Outgroup Interaction",
    y = "Outgroup Attitude",
    title = "Outgroup Attitudes by Interaction Type [95% CI]"
  ) +
  scale_fill_grey(
    start = 0.2,
    end = 0.8
  ) +
  scale_y_continuous(
    limits = c(0, 100),
    breaks = c(0, 15, 30, 40, 50, 60, 70, 85, 100),
    labels = c(
      "Very cold or unfavorable feelings 0°",
      "Quite cold and unfavorable feelings 15°",
      "Fairly cold and unfavorable feelings 30°",
      "A bit cold and unfavorable feelings 40°",
      "No feeling at all 50°",
      "A bit warm and favorable feelings 60°",
      "Fairly warm and favorable feelings 70° ",
      "Quite warm and favorable feelings 85° ",
      "Very warm and favorable feelings 100° "
    )
  ) +
  theme_Publication()

# create list to store Worker models
mdlWorker <- list()

# regression
mdlWorker$lmAttInt <-
  lm(thermometerDutch_1 ~ OutgroupInteraction * NonOutgroupInteraction,
    data = workerInteractionType
  )
# summary(lmWorkerAttInteraction)

summ(
  mdlWorker$lmAttInt,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type OLS linear regression
F(3,1221) 4.867
R² 0.012
Adj. R² 0.009
Est. 2.5% 97.5% t val. p
(Intercept) 69.507 67.899 71.116 84.796 0.000
OutgroupInteraction 3.620 0.378 6.862 2.191 0.029
NonOutgroupInteraction -0.513 -2.593 1.566 -0.484 0.628
OutgroupInteraction:NonOutgroupInteraction -0.098 -4.021 3.826 -0.049 0.961
Standard errors: OLS; Continuous predictors are mean-centered.

TL;DR: While controlling for interactions with non-Dutch people, outgroup attitudes were significantly higher when participants had an interaction with a Dutch person. The effect is relatively small (3.62 points on a 0–100 scale). More importantly, however, this analysis lumps all ESM beeps from every participants together and ignores that the data is nested within participants.

ML Regression with interaction term

Unconditional model

make a model with only the data but a random intercept (unconditional model)

\[\begin{equation} \tag{3} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmerAttNull <-
  lme4::lmer(thermometerDutch_1 ~ 1 + (1 | PID),
    data = dtWorker$full
  ) # use optim if it does not converge

mdlWorker$lmeAttNull <-
  lme(
    thermometerDutch_1 ~ 1,
    random = ~ 1 | PID,
    data = dtWorker$full,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(lmerWorkerAttNull) #or with the lme function
summ(mdlWorker$lmerAttNull, digits = 3, center = TRUE)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8805.880
BIC 8821.213
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.698
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 71.338 2.695 26.466 22.053 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.797
Residual 8.425
Grouping Variables
Group # groups ICC
PID 23 0.698 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorker$varAttNull <- 
  VarCorr(mdlWorker$lmeAttNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorker$tauAttNull <- 
  as.numeric(mdlWorker$varAttNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorker$sigmaAttNull <- 
  as.numeric(mdlWorker$varAttNull[2])
# The ICC estimate (between/between+within) is:
mdlWorker$IccAttNull <-
  (as.numeric(mdlWorker$varAttNull[1]) / (as.numeric(mdlWorker$varAttNull[1]) + as.numeric(mdlWorker$varAttNull[2])))
mdlWorker$IccPercAttNull <-
  ((as.numeric(mdlWorker$varAttNull[1]) / (as.numeric(mdlWorker$varAttNull[1]) + as.numeric(mdlWorker$varAttNull[2])))) * 100

ICC \(\tau_{00}^2 / (\tau_{00}^2 + \sigma^2)\): The ratio of the between-cluster variance to the total variance is called the Intraclass Correlation. It tells you the proportion of the total variance in Y that is accounted for by the clustering. (In this case the clustering means clustering observations per participant)

compare the random intercept model to a model without a random intercept (i.e., without levels at all)

# Create and save Model
mdlWorker$glsAttNull  <-
  gls(thermometerDutch_1 ~ 1,
      data = dtWorker$full,
      control = list(opt = "nlmimb"))

# calculate Deviances manually:
mdlWorker$DevianceGlsNull <- logLik(mdlWorker$glsAttNull) * -2
mdlWorker$DevianceMlNull <- logLik(mdlWorker$lmeAttNull) * -2

# Compare the two null models:
anova(mdlWorker$glsAttNull,
      mdlWorker$lmeAttNull) %>% 
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 1: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorker\(glsAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 2 </td> <td style="text-align:center;"> 10133 </td> <td style="text-align:center;"> 10144 </td> <td style="text-align:center;"> -5065 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlWorker\)lmeAttNull 2 3 8806 8821 -4400 1 vs 2 1329 0

We find that an estimated 69.76% of the variation in Feeling Thermometer scores is explained by between participant differences (clustering by PID). This is to say that 69.76% of the variance in any individual report of Attitudes towards the Dutch can be explained by the properties of the individual who provided the rating. And we find that including ‘participant’ as a predictor adds significantly to the model.

random intercept

\[\begin{equation} \tag{4} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmeInterceptAttType <-
  lme(
    thermometerDutch_1 ~ OutgroupInteraction + NonOutgroupInteraction,
    random =  ~ 1 | PID,
    data = workerInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    thermometerDutch_1 ~ OutgroupInteraction + NonOutgroupInteraction + (1 | PID),
    data = workerInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8788.218
BIC 8813.771
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.703
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.343 65.002 75.683 25.814 22.897 0.000
OutgroupInteraction 2.477 1.364 3.589 4.365 1204.135 0.000
NonOutgroupInteraction 0.427 -0.683 1.538 0.754 1204.911 0.451
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.811
Residual 8.362
Grouping Variables
Group # groups ICC
PID 23 0.701 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 2: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorker\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 8806 </td> <td style="text-align:center;"> 8821 </td> <td style="text-align:center;"> -4400 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlWorker\)lmeInterceptAttType 2 5 8788 8814 -4389 1 vs 2 21.66 0 
# Save variances
mdlWorker$varInterceptAttType <- 
  lme4::VarCorr(mdlWorker$lmeInterceptAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwInterceptAttType <-
  1 - (as.numeric(mdlWorker$varInterceptAttType[1]) / as.numeric(mdlWorker$varAttNull[1]))
mdlWorker$varBtwPercInterceptAttType <-
  (1 - (as.numeric(mdlWorker$varInterceptAttType[1]) / as.numeric(mdlWorker$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinInterceptAttType <-
  1 - (as.numeric(mdlWorker$varInterceptAttType[2]) / as.numeric(mdlWorker$varAttNull[2]))
mdlWorker$varWithinPercInterceptAttType <-
  (1 - (as.numeric(mdlWorker$varInterceptAttType[2]) / as.numeric(mdlWorker$varAttNull[2]))) * 100

tl;dr: Interaction with Dutch is great predictor, interactions with non-Dutch is not.

random slope

\[\begin{equation} \tag{5} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorker$lmeSlopesAttType <- lme(
  thermometerDutch_1 ~
    OutgroupInteraction + NonOutgroupInteraction,
  random = ~ 1 + OutgroupInteraction + NonOutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = workerInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorker$lmerSlopesAttType <- lmer(
    thermometerDutch_1 ~
      OutgroupInteraction + NonOutgroupInteraction +
      (1 + OutgroupInteraction + NonOutgroupInteraction | PID),
    data = workerInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8793.695
BIC 8844.802
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.710
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.371 65.001 75.741 25.682 21.938 0.000
OutgroupInteraction 2.572 1.079 4.065 3.376 19.291 0.003
NonOutgroupInteraction 0.434 -0.693 1.562 0.755 205.804 0.451
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.883
PID OutgroupInteraction 2.248
PID NonOutgroupInteraction 0.424
Residual 8.307
Grouping Variables
Group # groups ICC
PID 23 0.706 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeInterceptAttType, 
      mdlWorker$lmeSlopesAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 3: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorker\(lmeInterceptAttType </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 5 </td> <td style="text-align:center;"> 8788 </td> <td style="text-align:center;"> 8814 </td> <td style="text-align:center;"> -4389 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlWorker\)lmeSlopesAttType 2 10 8794 8845 -4387 1 vs 2 4.165 0.526 
# Save variances
mdlWorker$varSlopesAttType <- 
  lme4::VarCorr(mdlWorker$lmeSlopesAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwSlopesAttType <- 
  1 - (as.numeric(mdlWorker$varSlopesAttType[1]) / as.numeric(mdlWorker$varInterceptAttType[1]))
mdlWorker$varBtwPercSlopesAttType <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttType[1]) / as.numeric(mdlWorker$varInterceptAttType[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinSlopesAttType <- 
  1 - (as.numeric(mdlWorker$varSlopesAttType[2]) / as.numeric(mdlWorker$varInterceptAttType[2]))
mdlWorker$varWithinPercSlopesAttType <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttType[2]) / as.numeric(mdlWorker$varInterceptAttType[2]))) * 100

# Assumption Checks:
mdlWorker$diagSlopesAttType <-
  sjPlot::plot_model(mdlWorker$lmerSlopesAttType, type = "diag")
grid.arrange(
  mdlWorker$diagSlopesAttType[[1]],
  mdlWorker$diagSlopesAttType[[2]]$`PID`,
  mdlWorker$diagSlopesAttType[[3]],
  mdlWorker$diagSlopesAttType[[4]]
)

# Plot prediction model
mdlWorker$predSlopesAttType <- 
  workerInteractionType %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorker$lmeSlopesAttType,
                           workerInteractionType,
                           re.form = NA
                           )
         )

(
  mdlWorker$predPltSlopesAttType <-
    ggplot(data = mdlWorker$predSlopesAttType, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Worker_PredictionPlot_SlopesAttType.png",
  mdlWorker$predPltSlopesAttType,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

TL;DR: Random slopes don’t add much for this super simple model.

Interaction Frequency and Interaction Quality

\[\begin{equation} \tag{6} Attitude = ContactFreq \times AverageContactQual \end{equation}\]

# prepare data
workerAvFreQual <- dtWorker$full %>%
  group_by(ExternalReference) %>%
  summarise(
    n = n(),
    SumContactNL = sum(Contact_dum),
    PercContactNL = SumContactNL / n * 100,
    SumContactNLAll = sum(number),
    AvAttitude = mean(thermometerDutch_1, na.rm = TRUE),
    AvQuality = mean(quality_overall_1, na.rm = TRUE)
  ) %>%
  mutate(
    WinSumContactNL = DescTools::Winsorize(SumContactNL),
    WinSumContactNLAll = DescTools::Winsorize(SumContactNLAll)
  )

# correlation panel
pairs.panels.new(
  workerAvFreQual %>% select(SumContactNL, SumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  workerAvFreQual %>% select(WinSumContactNL, WinSumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

TL;DR: There is a medium sized correlation between the participants’ Average Interaction Quality and their Average Outgroup Attitudes. Thus within our data participants with a higher quality outgroup interactions also held more positive attitudes towards that group. However, the frequency of intergroup interactions had no meaningfull correlation with either the average interaction quality or average outgroup attitudes.

# regression
mdlWorker$lmAttFreqQualX <-
  lm(AvAttitude ~ SumContactNL * AvQuality, data = workerAvFreQual)
# summary(lmWorkerAttFreqQualX)

summ(
  mdlWorker$lmAttFreqQualX,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 21 (2 missing obs. deleted)
Dependent variable AvAttitude
Type OLS linear regression
F(3,17) 6.663
R² 0.540
Adj. R² 0.459
Est. 2.5% 97.5% t val. p
(Intercept) 70.930 66.519 75.341 33.927 0.000
SumContactNL 0.269 -0.150 0.688 1.354 0.193
AvQuality 0.765 0.288 1.242 3.385 0.004
SumContactNL:AvQuality -0.049 -0.085 -0.014 -2.954 0.009
Standard errors: OLS; Continuous predictors are mean-centered. 
interactions::interact_plot(
  mdlWorker$lmAttFreqQualX,
  pred = AvQuality,
  modx = SumContactNL,
  interval = TRUE,
  partial.residuals = TRUE
)

interactions::johnson_neyman(mdlWorker$lmAttFreqQualX,
                             pred = AvQuality,
                             modx = SumContactNL,
                             alpha = .05)
## JOHNSON-NEYMAN INTERVAL 
## 
## When SumContactNL is OUTSIDE the interval [23.66, 75.42], the slope of AvQuality is p < .05.
## 
## Note: The range of observed values of SumContactNL is [2.00, 51.00]

TL;DR: We find that interaction quality is significantly related to higher outgroup attitudes (albeit with a small effect size). We also find that in our sample with an increasing number of interactions the positive effect of interaction quality becomes weaker. However, note that this is based on data aggregating all within participant nuances and is only the date of 21 people.

Need Fulfillment

Need fulfillment and Attitudes

Unconditional model

make a model with only the data but a random intercept (unconditional model) for the outgroup interaction subset.

# see how large our outgroup interaction subset actually is
tbl_cross(
  workerInteractionType,
  row = OutgroupInteraction,
  col = NonOutgroupInteraction,
  percent = "row"
)
Characteristic Did you meet non-Dutch people this morning? (in person for at least 10 minutes) Total
no yes
Did you meet a Dutch person this morning? (In person interaction for at least 10 minutes)
No 337 (40%) 501 (60%) 838 (100%)
Yes 110 (28%) 277 (72%) 387 (100%)
Total 447 (36%) 778 (64%) 1,225 (100%)
 
# create outgroup interaction subset
workerOutgroupInteraction <- workerInteractionType %>%
  filter(OutgroupInteraction == "Yes")

# create empty list to organize models
mdlWorkerOut <- 
  list(
    Att = list(),
    Qlt = list()
  )

\[\begin{equation} \tag{7} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmerNull <-
  lme4::lmer(thermometerDutch_1 ~ 1 + (1 | PID), 
             data = workerOutgroupInteraction) # use optim if it does not converge
mdlWorkerOut$Att$lmeNull <-
  lme(
    thermometerDutch_1 ~ 1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.ML.r) #or with the lme function
summ(mdlWorkerOut$Att$lmerNull, digits = 3, center = TRUE)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2863.460
BIC 2875.336
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.684
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 72.550 2.910 24.933 19.198 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.069
Residual 8.883
Grouping Variables
Group # groups ICC
PID 21 0.684 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorkerOut$Att$varNull <- 
  VarCorr(mdlWorkerOut$Att$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorkerOut$Att$tauNull <- 
  as.numeric(mdlWorkerOut$Att$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorkerOut$Att$sigmaNull <- 
  as.numeric(mdlWorkerOut$Att$varNull[2])
# The ICC estimate (between/between+within) is:
mdlWorkerOut$Att$IccNull <-
  (as.numeric(mdlWorkerOut$Att$varNull[1]) / (as.numeric(mdlWorkerOut$Att$varNull[1]) + as.numeric(mdlWorkerOut$Att$varNull[2])))
mdlWorkerOut$Att$IccPercNull <-
  ((as.numeric(mdlWorkerOut$Att$varNull[1]) / (as.numeric(mdlWorkerOut$Att$varNull[1]) + as.numeric(mdlWorkerOut$Att$varNull[2])))) * 100

random intercept

\[\begin{equation} \tag{8} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCore <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(thermometerDutch_1 ~ keymotive_fulfillemt_1 + (1 | PID), 
       data = workerOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2842.675
BIC 2858.508
Pseudo-R² (fixed effects) 0.033
Pseudo-R² (total) 0.713
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.690 66.947 78.433 24.808 19.259 0.000
keymotive_fulfillemt_1 0.147 0.094 0.201 5.404 372.954 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.182
Residual 8.556
Grouping Variables
Group # groups ICC
PID 21 0.704 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeNull, 
      mdlWorkerOut$Att$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 4: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429 NA NA
mdlWorkerOut\(Att\)lmeInterceptCore 2 4 2843 2858 -1417 1 vs 2 22.79 0 
# Save variances
mdlWorkerOut$Att$varInterceptCore <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCore <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))
mdlWorkerOut$Att$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))
mdlWorkerOut$Att$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))) * 100

TL;DR: Situational core need relates significantly to outgroup attitudes. The core need has little explained variance between participants (compared to the null model: Variance Explained = 1 – (Var with Predictor/Var without Predictor); -1.73%). The variance explained within participants is small to medium (7.21%)

random slope

\[\begin{equation} \tag{9} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCore <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1,
    random = ~ 1 + keymotive_fulfillemt_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCore <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 +
      (1 + keymotive_fulfillemt_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2815.398
BIC 2839.148
Pseudo-R² (fixed effects) 0.046
Pseudo-R² (total) 0.763
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 71.863 66.133 77.594 24.578 19.326 0.000
keymotive_fulfillemt_1 0.179 0.061 0.297 2.970 18.186 0.008
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.111
PID keymotive_fulfillemt_1 0.219
Residual 7.941
Grouping Variables
Group # groups ICC
PID 21 0.732 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeInterceptCore, 
      mdlWorkerOut$Att$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 5: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeInterceptCore 1 4 2843 2858 -1417 NA NA
mdlWorkerOut\(Att\)lmeSlopesCore 2 6 2815 2839 -1402 1 vs 2 31.28 0 
# Save variances
mdlWorkerOut$Att$varSlopesCore <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))
mdlWorkerOut$Att$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))
mdlWorkerOut$Att$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCore <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCore[[1]],
  mdlWorkerOut$Att$diagSlopesCore[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCore[[3]],
  mdlWorkerOut$Att$diagSlopesCore[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCore <- 
  workerOutgroupInteraction %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCore,
                           workerOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCore <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCore, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCore.png",
  mdlWorkerOut$Att$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

TL;DR: The random slope adds significantly to the prediction model.

Need fulfillment and Interaction Quality

Unconditional model

make a model with only the data but a random intercept (unconditional model) for the outgroup interaction subset.

\[\begin{equation} \tag{10} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Qlt$lmerNull <-
  lme4::lmer(quality_overall_1 ~ 1 + (1 | PID), 
             data = workerOutgroupInteraction) # use optim if it does not converge
mdlWorkerOut$Qlt$lmeNull <-
  lme(
    quality_overall_1 ~ 1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.Qual.ML.r) #or with the lme function
summ(mdlWorkerOut$Qlt$lmerNull, digits = 3, center = TRUE)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3347.534
BIC 3359.410
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.183
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 24.285 2.090 11.619 20.156 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 8.286
Residual 17.511
Grouping Variables
Group # groups ICC
PID 21 0.183 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1~1 + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlWorkerOut$Qlt$varNull <- 
  VarCorr(mdlWorkerOut$Qlt$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlWorkerOut$Qlt$tauNull <- 
  as.numeric(mdlWorkerOut$Qlt$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlWorkerOut$Qlt$sigmaNull <- 
  as.numeric(mdlWorkerOut$Qlt$varNull[2])
# The ICC estimate (between/between+within) is:
mdlWorkerOut$Qlt$IccNull <-
  (as.numeric(mdlWorkerOut$Qlt$varNull[1]) / (as.numeric(mdlWorkerOut$Qlt$varNull[1]) + as.numeric(mdlWorkerOut$Qlt$varNull[2])))
mdlWorkerOut$Qlt$IccPercNull <-
  ((as.numeric(mdlWorkerOut$Qlt$varNull[1]) / (as.numeric(mdlWorkerOut$Qlt$varNull[1]) + as.numeric(mdlWorkerOut$Qlt$varNull[2])))) * 100

random intercept

\[\begin{equation} \tag{11} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Qlt$lmeInterceptCore <-
  lme(
    quality_overall_1 ~ keymotive_fulfillemt_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(quality_overall_1 ~ keymotive_fulfillemt_1 + (1 | PID), 
       data = workerOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3287.252
BIC 3303.085
Pseudo-R² (fixed effects) 0.179
Pseudo-R² (total) 0.303
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 24.551 21.086 28.016 13.886 19.850 0.000
keymotive_fulfillemt_1 0.417 0.321 0.513 8.528 360.783 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.814
Residual 16.157
Grouping Variables
Group # groups ICC
PID 21 0.151 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(quality_overall_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Qlt$lmeNull, 
      mdlWorkerOut$Qlt$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 6: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Qlt\)lmeNull 1 3 3348 3359 -1671 NA NA
mdlWorkerOut\(Qlt\)lmeInterceptCore 2 4 3287 3303 -1640 1 vs 2 62.28 0 
# Save variances
mdlWorkerOut$Qlt$varInterceptCore <-
  lme4::VarCorr(mdlWorkerOut$Qlt$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Qlt$varBtwInterceptCore <- 
  1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Qlt$varNull[1]))
mdlWorkerOut$Qlt$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]) / as.numeric(mdlWorkerOut$Qlt$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Qlt$varWithinInterceptCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Qlt$varNull[2]))
mdlWorkerOut$Qlt$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]) / as.numeric(mdlWorkerOut$Qlt$varNull[2]))) * 100

TL;DR: Situational core need relates significantly to interaction quality. The core need explained substantial variance between participants (32.35%). The variance explained within participants is also medium (14.87%)

random slope

\[\begin{equation} \tag{12} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Qlt$lmeSlopesCore <-
  lme(
    quality_overall_1 ~
      keymotive_fulfillemt_1,
    random = ~ 1 + keymotive_fulfillemt_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Qlt$lmerSlopesCore <-
    lmer(
      quality_overall_1 ~
        keymotive_fulfillemt_1 +
        (1 + keymotive_fulfillemt_1 | PID),
      data = workerOutgroupInteraction
    ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable quality_overall_1
Type Mixed effects linear regression
AIC 3289.190
BIC 3312.941
Pseudo-R² (fixed effects) 0.192
Pseudo-R² (total) 0.357
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 24.263 20.622 27.903 13.063 18.950 0.000
keymotive_fulfillemt_1 0.443 0.304 0.582 6.250 7.336 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 7.176
PID keymotive_fulfillemt_1 0.189
Residual 15.918
Grouping Variables
Group # groups ICC
PID 21 0.169 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(mdlWorkerOut$Qlt$lmerSlopesCore,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Qlt$lmeInterceptCore, 
      mdlWorkerOut$Qlt$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 7: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Qlt\)lmeInterceptCore 1 4 3287 3303 -1640 NA NA
mdlWorkerOut\(Qlt\)lmeSlopesCore 2 6 3289 3313 -1639 1 vs 2 2.062 0.357 
# Save variances
mdlWorkerOut$Qlt$varSlopesCore <- 
  lme4::VarCorr(mdlWorkerOut$Qlt$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Qlt$varBtwSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]))
mdlWorkerOut$Qlt$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[1]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Qlt$varWithinSlopesCore <-
  1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]))
mdlWorkerOut$Qlt$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlWorkerOut$Qlt$varSlopesCore[2]) / as.numeric(mdlWorkerOut$Qlt$varInterceptCore[2]))) * 100

TL;DR: random slopes does not add significantly. Stay with simpler model.

Need fulfillment, Interaction Quality, and Attitudes

random intercept

\[\begin{equation} \tag{13} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCoreQlt <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + quality_overall_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + quality_overall_1 + (1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2756.572
BIC 2776.364
Pseudo-R² (fixed effects) 0.123
Pseudo-R² (total) 0.751
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.697 67.464 77.931 27.225 19.335 0.000
keymotive_fulfillemt_1 0.042 -0.009 0.094 1.614 371.091 0.107
quality_overall_1 0.252 0.204 0.300 10.288 367.215 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.025
Residual 7.574
Grouping Variables
Group # groups ICC
PID 21 0.716 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(
  mdlWorkerOut$Att$lmeNull, 
  mdlWorkerOut$Att$lmeInterceptCoreQlt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 8: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429 NA NA
mdlWorkerOut\(Att\)lmeInterceptCoreQlt 2 5 2757 2776 -1373 1 vs 2 110.9 0 
# Save variances
mdlWorkerOut$Att$varInterceptCoreQlt <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCoreQlt <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))
mdlWorkerOut$Att$varBtwPercInterceptCoreQlt <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))
mdlWorkerOut$Att$varWithinPercInterceptCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varNull[2]))) * 100

TL;DR: Variance explained in outgroup attitudes is of medium effect size (between: 15.34%, within: 27.29%)

random slope

\[\begin{equation} \tag{14} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCoreQlt <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + quality_overall_1,
    random = ~ 1 + keymotive_fulfillemt_1 + quality_overall_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCoreQlt <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + quality_overall_1 +
      (1 + keymotive_fulfillemt_1 + quality_overall_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 387
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2669.434
BIC 2709.018
Pseudo-R² (fixed effects) 0.105
Pseudo-R² (total) 0.836
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 72.553 67.394 77.711 27.566 19.078 0.000
keymotive_fulfillemt_1 0.036 -0.076 0.149 0.636 15.027 0.535
quality_overall_1 0.233 0.128 0.338 4.345 21.118 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 11.827
PID keymotive_fulfillemt_1 0.216
PID quality_overall_1 0.210
Residual 6.133
Grouping Variables
Group # groups ICC
PID 21 0.788 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(
  mdlWorkerOut$Att$lmeNull,
  mdlWorkerOut$Att$lmeInterceptCoreQlt,
  mdlWorkerOut$Att$lmeSlopesCoreQlt
) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 9: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeNull 1 3 2863 2875 -1429 NA NA
mdlWorkerOut\(Att\)lmeInterceptCoreQlt 2 5 2757 2776 -1373 1 vs 2 110.89 0
mdlWorkerOut\(Att\)lmeSlopesCoreQlt 3 10 2669 2709 -1325 2 vs 3 97.14 0 
# Save variances
mdlWorkerOut$Att$varSlopesCoreQlt <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]))
mdlWorkerOut$Att$varBtwPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCoreQlt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]))
mdlWorkerOut$Att$varWithinPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreQlt[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCoreQlt <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCoreQlt, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCoreQlt[[1]],
  mdlWorkerOut$Att$diagSlopesCoreQlt[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCoreQlt[[3]],
  mdlWorkerOut$Att$diagSlopesCoreQlt[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCoreQlt <- 
  workerOutgroupInteraction %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCoreQlt,
                           workerOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCoreQlt <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCoreQlt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCoreQlt.png",
  mdlWorkerOut$Att$predPltSlopesCoreQlt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

Compare models

Interpretation of Mediation like results.

Check for robustness

Interaction Type

random intercept

\[\begin{equation} \tag{15} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorker$lmeInterceptAttCoreInt <-
  lme(
    thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction,
    random =  ~ 1 | PID,
    data = workerInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlWorker$lmerInterceptAttCoreInt <- lmer(
    thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction + (1 | PID),
    data = workerInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8774.952
BIC 8805.616
Pseudo-R² (fixed effects) 0.015
Pseudo-R² (total) 0.713
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.587 65.258 75.916 25.961 22.187 0.000
keyMotiveFulfilled 0.016 -0.007 0.039 1.377 1205.341 0.169
OutgroupInteraction 1.711 0.578 2.843 2.961 1202.807 0.003
keyMotiveFulfilled:OutgroupInteraction 0.109 0.061 0.157 4.428 1200.535 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.891
Residual 8.265
Grouping Variables
Group # groups ICC
PID 23 0.709 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=dtWorker$full),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 10: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorker\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 8806 </td> <td style="text-align:center;"> 8821 </td> <td style="text-align:center;"> -4400 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlWorker\)lmeInterceptAttCoreInt 2 6 8775 8806 -4381 1 vs 2 36.93 0 
# Save variances
mdlWorker$varInterceptAttCoreInt <- 
  lme4::VarCorr(mdlWorker$lmeInterceptAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwInterceptAttCoreInt <-
  1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[1]) / as.numeric(mdlWorker$varAttNull[1]))
mdlWorker$varBtwPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[1]) / as.numeric(mdlWorker$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinInterceptAttCoreInt <-
  1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[2]) / as.numeric(mdlWorker$varAttNull[2]))
mdlWorker$varWithinPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlWorker$varInterceptAttCoreInt[2]) / as.numeric(mdlWorker$varAttNull[2]))) * 100
random slope

\[\begin{equation} \tag{16} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorker$lmeSlopesAttCoreInt <- lme(
  thermometerDutch_1 ~
    keyMotiveFulfilled * OutgroupInteraction,
  random = ~ 1 + keyMotiveFulfilled + OutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = workerInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorker$lmerSlopesAttCoreInt <- lmer(
    thermometerDutch_1 ~
      keyMotiveFulfilled * OutgroupInteraction +
      (1 + keyMotiveFulfilled + OutgroupInteraction | PID),
    data = workerInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 1225
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 8763.482
BIC 8819.700
Pseudo-R² (fixed effects) 0.020
Pseudo-R² (total) 0.729
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.252 65.044 75.459 26.440 22.117 0.000
keyMotiveFulfilled 0.031 -0.021 0.082 1.166 10.093 0.270
OutgroupInteraction 1.846 0.169 3.523 2.157 19.514 0.044
keyMotiveFulfilled:OutgroupInteraction 0.112 0.057 0.167 3.984 396.330 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 12.558
PID keyMotiveFulfilled 0.100
PID OutgroupInteraction 2.795
Residual 8.043
Grouping Variables
Group # groups ICC
PID 23 0.709 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorker$lmeAttNull, 
      mdlWorker$lmeInterceptAttCoreInt,
      mdlWorker$lmeSlopesAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 11: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorker\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 8806 </td> <td style="text-align:center;"> 8821 </td> <td style="text-align:center;"> -4400 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlWorker\)lmeInterceptAttCoreInt 2 6 8775 8806 -4381 1 vs 2 36.93 0.000
mdlWorker$lmeSlopesAttCoreInt 3 11 8763 8820 -4371 2 vs 3 21.47 0.001 
# Save variances
mdlWorker$varSlopesAttCoreInt <- 
  lme4::VarCorr(mdlWorker$lmeSlopesAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorker$varBtwSlopesAttCoreInt <- 
  1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[1]) / as.numeric(mdlWorker$varInterceptAttCoreInt[1]))
mdlWorker$varBtwPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[1]) / as.numeric(mdlWorker$varInterceptAttCoreInt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorker$varWithinSlopesAttCoreInt <- 
  1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[2]) / as.numeric(mdlWorker$varInterceptAttCoreInt[2]))
mdlWorker$varWithinPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlWorker$varSlopesAttCoreInt[2]) / as.numeric(mdlWorker$varInterceptAttCoreInt[2]))) * 100

# Assumption Checks:
mdlWorker$diagSlopesAttCoreInt <-
  sjPlot::plot_model(mdlWorker$lmerSlopesAttCoreInt, type = "diag")
grid.arrange(
  mdlWorker$diagSlopesAttCoreInt[[1]],
  mdlWorker$diagSlopesAttCoreInt[[2]]$`PID`,
  mdlWorker$diagSlopesAttCoreInt[[3]],
  mdlWorker$diagSlopesAttCoreInt[[4]]
)

# Plot prediction model
mdlWorker$predSlopesAttCoreInt <- 
  workerInteractionType %>%
  select(thermometerDutch_1, session, PID) %>% 
  mutate(measure = predict(mdlWorker$lmeSlopesAttCoreInt,
                           workerInteractionType,
                           re.form = NA
                           )
         )

(
  mdlWorker$predPltSlopesAttCoreInt <-
    ggplot(data = mdlWorker$predSlopesAttCoreInt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Worker_PredictionPlot_SlopesAttCoreInt.png",
  mdlWorker$predPltSlopesAttCoreInt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)
Plot Interaction
# visualize interaction
## Without ML structure
ggplot(data = workerInteractionType,
       aes(x = keyMotiveFulfilled,
           y = thermometerDutch_1,
           fill = OutgroupInteraction)) +
  #geom_point()+
  geom_smooth(method = 'lm',
              aes(linetype = OutgroupInteraction),
              color = "black") +
  #facet_wrap(~PID, ncol = 6)+
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  #scale_colour_manual(values=c("grey20", "black"), name="Intergroup Contact")+
  scale_y_continuous(
    limits = c(50, 100),
    breaks = seq(50, 100, by = 10),
    position = "left"
  ) +
  scale_x_continuous(limits = c(-50, 50), breaks = seq(-50, 50, by = 10)) +
  labs(
    title = "Without ML stucture",
    x = "Fulfillment Core Need",
    y = "Outgroup Attitudes",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact"
  ) +
  theme_Publication() +
  theme(legend.position = "bottom",
        legend.key.size = unit(1, "cm"))

## With ML structure
# create parameters for prediction
datNew = data.frame(
  OutgroupInteraction = factor(rep(c("No", "Yes"), each = 21)),
  keyMotiveFulfilled = rep(seq(-50, 50, 5), 2),
  PID = 0
)

# Predict values, clean up and calculate SE
PI <-
  merTools::predictInterval(
    merMod = mdlWorker$lmerSlopesAttCoreInt,
    newdata = datNew,
    level = 0.95,
    stat = "mean",
    type = "linear.prediction",
    include.resid.var = F,
    fix.intercept.variance = F
  )
mdlWorker$predInterceptAttCoreIntX <- 
  cbind(datNew, PI)
mdlWorker$predInterceptAttCoreIntX$se <-
  (mdlWorker$predInterceptAttCoreIntX$upr - mdlWorker$predInterceptAttCoreIntX$fit) / 1.96
rm(datNew, PI)

# Plot predicted values with SE
ggplot(
  mdlWorker$predInterceptAttCoreIntX,
  aes(x = keyMotiveFulfilled,
      y = fit,
      fill = OutgroupInteraction)
)+
  #geom_point() +
  geom_line(aes(linetype = as.factor(OutgroupInteraction)), size = 1) +
  #facet_wrap(~PID, ncol = 6)+
  geom_ribbon(data = mdlWorker$predInterceptAttCoreIntX,
              aes(ymin = fit - se, ymax = fit + se),
              alpha = 0.3) +
  scale_x_continuous(breaks = seq(-50, 50, 10)) +
  scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  labs(
    x = "Fulfillment Core Need",
    y = "Outgroup Attitude (NL)",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact",
    title = "Based on Model Predictions"
  ) +
  theme_Publication()

# #### Bayesian estimation !! ONLY RUN ON FINAL RENDER !! Takes forever ####
# options(mc.cores = parallel::detectCores())  # Run many chains simultaneously
# brmfit <- brm(
#   thermometerDutch_1 ~ keyMotiveFulfilled * OutgroupInteraction + 
#     (1 + keyMotiveFulfilled + OutgroupInteraction | PID),
#   data = workerInteractionType,
#   family = gaussian,
#   iter = 1000,
#   chains = 4
# )
# 
# # create parameters for prediction:
# datNew = data.frame(
#   OutgroupInteraction = rep(c("No", "Yes"), each = 101),
#   keyMotiveFulfilled = rep(seq(-50, 50, 1), 2)
# )
# # Save predicted values and adjust names and labels
# fitavg <-
#   cbind(datNew,
#         fitted(brmfit, newdata = datNew, re_formula = NA)[, -2])
# names(fitavg)[names(fitavg) == "Estimate"] = "pred"
# fitavg$se <- (fitavg$Q97.5 - fitavg$pred) / 1.96
# 
# # Plot Bayesian SE prediction interval
# ggplot(fitavg,
#        aes(x = keyMotiveFulfilled,
#            y = pred,
#            fill = OutgroupInteraction)) +
#   scale_x_continuous(breaks = seq(-50, 50, 10)) +
#   scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
#   geom_line(aes(linetype = OutgroupInteraction), size = 1) +
#   geom_ribbon(aes(ymin = pred - se, ymax = pred + se), alpha = 0.2) +
#   scale_linetype_manual(values = c("dashed", "solid")) +
#   scale_fill_manual(values = c("darkgrey", "black")) +
#   labs(
#     x = "Fulfillment Core Need",
#     y = "Outgroup Attitude (NL)",
#     fill = "Intergroup Contact",
#     linetype = "Intergroup Contact",
#     title = "Based on Bayesian Prediction Interval"
#   ) +
#   theme_Publication()
# 
# # plot all overlayed posteriors:
# pst <- posterior_samples(brmfit, "b")
# ggplot(workerInteractionType,
#        aes(x = keyMotiveFulfilled, y = thermometerDutch_1)) +
#   geom_point(shape = 4, alpha = .1) +
#   geom_tile() +
#   geom_abline(
#     data = pst,
#     aes(intercept = b_Intercept, slope = b_keyMotiveFulfilled),
#     alpha = .025,
#     size = .4
#   ) +
#   labs(title = "slope Posteriors",
#        x = "Fulfillment Core Need",
#        y = "Outgroup Attitudes (NL)") +
#   theme_Publication()
# rm(datNew, brmfit, fitavg, pst)

Other psychological needs

random intercept

\[\begin{equation} \tag{17} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlWorkerOut$Att$lmeInterceptCoreSdt <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlWorkerOut$Att$lmerInterceptCoreSdt <- lmer(
    thermometerDutch_1 ~ keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 + (1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 386
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2804.147
BIC 2831.838
Pseudo-R² (fixed effects) 0.086
Pseudo-R² (total) 0.754
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 73.077 67.364 78.790 25.071 19.379 0.000
keymotive_fulfillemt_1 0.092 0.038 0.146 3.323 365.656 0.001
competence_1 0.019 -0.034 0.072 0.705 363.428 0.481
autonomy_1 0.079 0.014 0.143 2.388 363.976 0.017
relatedness_1 0.114 0.078 0.150 6.204 365.509 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.139
Residual 7.978
Grouping Variables
Group # groups ICC
PID 21 0.731 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(thermometerDutch_1 ~ Contact_dum + inNonDutch + (1|PID),data=workerOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# compute lme for core model with same number of observations for comparison
# !! LOOK INTO MISSINGNESS (only one observation) !!
lmeInterceptCoreRed <-
  lme(
    thermometerDutch_1 ~ keymotive_fulfillemt_1,
    random = ~ 1 | PID,
    data = workerOutgroupInteraction %>% filter(!is.na(autonomy_1))
  )

# Compare new model to previous step
anova(
  lmeInterceptCoreRed, 
  mdlWorkerOut$Att$lmeInterceptCoreSdt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 12: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
lmeInterceptCoreRed 1 4 2836 2852 -1414 NA NA
mdlWorkerOut\(Att\)lmeInterceptCoreSdt 2 7 2804 2832 -1395 1 vs 2 38.22 0 
rm(lmeInterceptCoreRed)

# Save variances
mdlWorkerOut$Att$varInterceptCoreSdt <-
  lme4::VarCorr(mdlWorkerOut$Att$lmeInterceptCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwInterceptCoreSdt <- 
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))
mdlWorkerOut$Att$varBtwPercInterceptCoreSdt <- 
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinInterceptCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))
mdlWorkerOut$Att$varWithinPercInterceptCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCore[2]))) * 100
random slope

\[\begin{equation} \tag{18} Attitude_{ti} = \\gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlWorkerOut$Att$lmeSlopesCoreSdt <-
  lme(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1,
    random = ~ 1 + keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 | PID,
    control = lmeControl(opt = "optim"),
    data = workerOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlWorkerOut$Att$lmerSlopesCoreSdt <- lmer(
    thermometerDutch_1 ~
      keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 +
      (1 + keymotive_fulfillemt_1 + competence_1 + autonomy_1 + relatedness_1 | PID),
    data = workerOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 386
Dependent variable thermometerDutch_1
Type Mixed effects linear regression
AIC 2754.434
BIC 2837.507
Pseudo-R² (fixed effects) 0.067
Pseudo-R² (total) 0.870
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.956 65.168 76.745 24.027 18.872 0.000
keymotive_fulfillemt_1 0.091 0.012 0.169 2.259 10.955 0.045
competence_1 0.047 -0.166 0.259 0.429 9.300 0.678
autonomy_1 0.059 -0.120 0.239 0.646 9.357 0.534
relatedness_1 0.100 0.056 0.145 4.412 18.715 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 13.289
PID keymotive_fulfillemt_1 0.123
PID competence_1 0.467
PID autonomy_1 0.378
PID relatedness_1 0.064
Residual 6.423
Grouping Variables
Group # groups ICC
PID 21 0.811 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlWorkerOut$Att$lmeInterceptCoreSdt, 
      mdlWorkerOut$Att$lmeSlopesCoreSdt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Worker: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 13: Worker: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlWorkerOut\(Att\)lmeInterceptCoreSdt 1 7 2804 2832 -1395 NA NA
mdlWorkerOut\(Att\)lmeSlopesCoreSdt 2 21 2756 2838 -1357 1 vs 2 76.51 0 
# Save variances
mdlWorkerOut$Att$varSlopesCoreSdt <- 
  lme4::VarCorr(mdlWorkerOut$Att$lmeSlopesCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlWorkerOut$Att$varBtwSlopesCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]))
mdlWorkerOut$Att$varBtwPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlWorkerOut$Att$varWithinSlopesCoreSdt <-
  1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]))
mdlWorkerOut$Att$varWithinPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlWorkerOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlWorkerOut$Att$varInterceptCoreSdt[2]))) * 100

# Assumption Checks:
mdlWorkerOut$Att$diagSlopesCoreSdt <- 
  sjPlot::plot_model(mdlWorkerOut$Att$lmerSlopesCoreSdt, type = "diag")
grid.arrange(
  mdlWorkerOut$Att$diagSlopesCoreSdt[[1]],
  mdlWorkerOut$Att$diagSlopesCoreSdt[[2]]$`PID`,
  mdlWorkerOut$Att$diagSlopesCoreSdt[[3]],
  mdlWorkerOut$Att$diagSlopesCoreSdt[[4]]
)

# Plot prediction model
mdlWorkerOut$Att$predSlopesCoreSdt <- 
  workerOutgroupInteraction %>%
  filter(!is.na(autonomy_1)) %>%
  select(thermometerDutch_1, session, PID, autonomy_1) %>% 
  mutate(measure = predict(mdlWorkerOut$Att$lmeSlopesCoreSdt,
                           workerOutgroupInteraction %>% filter(!is.na(autonomy_1)),
                           re.form = NA
                           )
         )

(
  mdlWorkerOut$Att$predPltSlopesCoreSdt <-
    ggplot(data = mdlWorkerOut$Att$predSlopesCoreSdt, aes(x = session, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = thermometerDutch_1), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/WorkerOut_PredictionPlot_SlopesAttCoreStd.png",
  mdlWorkerOut$Att$predPltSlopesCoreSdt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

Student Sample

Data Description

Still in ‘scr/studentDescriptive.R’. Needs to be merged with this document.

Contact Hypothesis

Interaction Frequency and Outgroup Attitudes

\[\begin{equation} \tag{19} r_{ContactFrequency, OutgroupAttitudes} \neq 0 \end{equation}\]

# correlation panel
pairs.panels.new(
  studentContactFreq %>% select(SumContactNL, SumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  studentContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

TL;DR: Neither the number of interactions nor the number of measurement beeps with an interaction are significantly related with the average outgroup attitudes. This is to say that within our data participants with more outgroup interactions did not have significantly more positive outgroup attitudes. This might be due to the aggregation within the participants or the small sample size of between participant data. Nonetheless, the aggregate data does not support the notion that simply having more interactions with an outgroup results in more positive outgroup attitudes.

Outgroup Attitudes by Interaction Type

\[\begin{equation} \tag{20} \mu_{i,0utgroupInteraction} > \mu_{i,IngroupInteraction} \\ Attitude = OutgroupInteraction \times NonOutgroupInteraction \end{equation}\]

Regression with interaction term

Across participants: impact of different interaction types [probably meaningless because ignores response sets between participants]

# between participants interaction type
studentContactType <- studentInteractionType %>%
  group_by(
    OutgroupInteraction,
    NonOutgroupInteraction
  ) %>%
  summarise(
    n = n(),
    AttitudeM = mean(AttitudesDutch, na.rm = TRUE),
    AttitudeSD = sd(AttitudesDutch, na.rm = TRUE),
    AttitudeSE = AttitudeSD / sqrt(n),
    AttitudeLwr = AttitudeM - 1.96 * AttitudeSE,
    AttitudeUpr = AttitudeM + 1.96 * AttitudeSE
  ) %>%
  ungroup()

# plot bar chart
ggplot(
  studentContactType,
  aes(
    y = AttitudeM,
    x = OutgroupInteraction,
    fill = NonOutgroupInteraction
  )
) +
  geom_bar(
    stat = "identity",
    color = "black",
    position = position_dodge()
  ) +
  geom_errorbar(aes(ymin = AttitudeM, ymax = AttitudeUpr),
    width = .2,
    position = position_dodge(.9)
  ) +
  labs(
    fill = "Non-Outgroup Interaction",
    x = "Outgroup Interaction",
    y = "Outgroup Attitude",
    title = "Outgroup Attitudes by Interaction Type [95% CI]"
  ) +
  scale_fill_grey(
    start = 0.2,
    end = 0.8
  ) +
  scale_y_continuous(
    limits = c(0, 100),
    breaks = c(0, 15, 30, 40, 50, 60, 70, 85, 100),
    labels = c(
      "Very cold or unfavorable feelings 0°",
      "Quite cold and unfavorable feelings 15°",
      "Fairly cold and unfavorable feelings 30°",
      "A bit cold and unfavorable feelings 40°",
      "No feeling at all 50°",
      "A bit warm and favorable feelings 60°",
      "Fairly warm and favorable feelings 70° ",
      "Quite warm and favorable feelings 85° ",
      "Very warm and favorable feelings 100° "
    )
  ) +
  theme_Publication()

# create list to store student models
mdlStudent <- list()

# regression
mdlStudent$lmAttInt <-
  lm(AttitudesDutch ~ OutgroupInteraction * NonOutgroupInteraction,
    data = studentInteractionType
  )
# summary(lmstudentAttInteraction)

summ(
  mdlStudent$lmAttInt,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type OLS linear regression
F(3,4961) 37.435
R² 0.022
Adj. R² 0.022
Est. 2.5% 97.5% t val. p
(Intercept) 66.807 65.815 67.799 132.053 0.000
OutgroupInteraction 8.159 5.699 10.619 6.503 0.000
NonOutgroupInteraction -1.226 -2.529 0.077 -1.845 0.065
OutgroupInteraction:NonOutgroupInteraction -0.355 -3.439 2.730 -0.225 0.822
Standard errors: OLS; Continuous predictors are mean-centered.

TL;DR: While controlling for interactions with non-Dutch people, outgroup attitudes were significantly higher when participants had an interaction with a Dutch person. The effect is relatively small (8.16 points on a 0–100 scale). More importantly, however, this analysis lumps all ESM beeps from every participants together and ignores that the data is nested within participants.

ML Regression with interaction term

Unconditional model

make a model with only the data but a random intercept (unconditional model)

\[\begin{equation} \tag{21} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmerAttNull <-
  lme4::lmer(AttitudesDutch ~ 1 + (1 | PID),
    data = dtStudents$full
  ) # use optim if it does not converge

mdlStudent$lmeAttNull <-
  lme(
    AttitudesDutch ~ 1,
    random = ~ 1 | PID,
    data = dtStudents$full,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(mdlStudent$lmerAttNull) #or with the lme function
summ(mdlStudent$lmerAttNull, digits = 3, center = TRUE)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36740.926
BIC 36760.456
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.801
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 67.290 1.751 38.435 111.704 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.540
Residual 9.235
Grouping Variables
Group # groups ICC
PID 113 0.801 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=dtStudent$full),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudent$varAttNull <- 
  VarCorr(mdlStudent$lmeAttNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudent$tauAttNull <- 
  as.numeric(mdlStudent$varAttNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudent$sigmaAttNull <- 
  as.numeric(mdlStudent$varAttNull[2])
# The ICC estimate (between/between+within) is:
mdlStudent$IccAttNull <-
  (as.numeric(mdlStudent$varAttNull[1]) / (as.numeric(mdlStudent$varAttNull[1]) + as.numeric(mdlStudent$varAttNull[2])))
mdlStudent$IccPercAttNull <-
  ((as.numeric(mdlStudent$varAttNull[1]) / (as.numeric(mdlStudent$varAttNull[1]) + as.numeric(mdlStudent$varAttNull[2])))) * 100

ICC \(\tau_{00}^2 / (\tau_{00}^2 + \sigma^2)\): The ratio of the between-cluster variance to the total variance is called the Intraclass Correlation. It tells you the proportion of the total variance in Y that is accounted for by the clustering. (In this case the clustering means clustering observations per participant)

compare the random intercept model to a model without a random intercept (i.e., without levels at all)

# Create and save Model
mdlStudent$glsAttNull  <-
  gls(AttitudesDutch ~ 1,
      data = dtStudents$full,
      control = list(opt = "nlmimb"))

# calculate Deviances manually:
mdlStudent$DevianceGlsNull <- logLik(mdlStudent$glsAttNull) * -2
mdlStudent$DevianceMlNull <- logLik(mdlStudent$lmeAttNull) * -2

# Compare the two null models:
anova(mdlStudent$glsAttNull,
      mdlStudent$lmeAttNull) %>% 
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 14: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudent\(glsAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 2 </td> <td style="text-align:center;"> 44352 </td> <td style="text-align:center;"> 44365 </td> <td style="text-align:center;"> -22174 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlStudent\)lmeAttNull 2 3 36741 36760 -18367 1 vs 2 7614 0

We find that an estimated 80.12% of the variation in Feeling Thermometer scores is explained by between participant differences (clustering by PID). This is to say that 80.12% of the variance in any individual report of Attitudes towards the Dutch can be explained by the properties of the individual who provided the rating. And we find that including ‘participant’ as a predictor adds significantly to the model.

random intercept

\[\begin{equation} \tag{22} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmeInterceptAttType <-
  lme(
    AttitudesDutch ~ OutgroupInteraction + NonOutgroupInteraction,
    random =  ~ 1 | PID,
    data = studentInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    AttitudesDutch ~ OutgroupInteraction + NonOutgroupInteraction + (1 | PID),
    data = studentInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36684.392
BIC 36716.943
Pseudo-R² (fixed effects) 0.003
Pseudo-R² (total) 0.801
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.800 63.374 70.226 38.217 114.319 0.000
OutgroupInteraction 2.882 2.162 3.601 7.852 4862.093 0.000
NonOutgroupInteraction -0.133 -0.708 0.443 -0.453 4861.600 0.651
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.400
Residual 9.180
Grouping Variables
Group # groups ICC
PID 113 0.801 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentInteractionType),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 15: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudent\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 36741 </td> <td style="text-align:center;"> 36760 </td> <td style="text-align:center;"> -18367 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlStudent\)lmeInterceptAttType 2 5 36684 36717 -18337 1 vs 2 60.53 0 
# Save variances
mdlStudent$varInterceptAttType <- 
  lme4::VarCorr(mdlStudent$lmeInterceptAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwInterceptAttType <-
  1 - (as.numeric(mdlStudent$varInterceptAttType[1]) / as.numeric(mdlStudent$varAttNull[1]))
mdlStudent$varBtwPercInterceptAttType <-
  (1 - (as.numeric(mdlStudent$varInterceptAttType[1]) / as.numeric(mdlStudent$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinInterceptAttType <-
  1 - (as.numeric(mdlStudent$varInterceptAttType[2]) / as.numeric(mdlStudent$varAttNull[2]))
mdlStudent$varWithinPercInterceptAttType <-
  (1 - (as.numeric(mdlStudent$varInterceptAttType[2]) / as.numeric(mdlStudent$varAttNull[2]))) * 100

tl;dr: Interaction with Dutch is great predictor, interactions with non-Dutch is not.

random slope

\[\begin{equation} \tag{5} Attitude_{ti} = \gamma_{00} + \gamma_{10}OutgroupInteraction_{ti} + \gamma_{10}NonOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudent$lmeSlopesAttType <- lme(
  AttitudesDutch ~
    OutgroupInteraction + NonOutgroupInteraction,
  random = ~ 1 + OutgroupInteraction + NonOutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = studentInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudent$lmerSlopesAttType <- lmer(
    AttitudesDutch ~
      OutgroupInteraction + NonOutgroupInteraction +
      (1 + OutgroupInteraction + NonOutgroupInteraction | PID),
    data = studentInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36476.495
BIC 36541.597
Pseudo-R² (fixed effects) 0.003
Pseudo-R² (total) 0.817
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.717 63.293 70.141 38.191 111.810 0.000
OutgroupInteraction 2.992 1.448 4.537 3.797 93.820 0.000
NonOutgroupInteraction 0.009 -0.763 0.781 0.022 108.800 0.982
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.375
PID OutgroupInteraction 6.996
PID NonOutgroupInteraction 2.689
Residual 8.801
Grouping Variables
Group # groups ICC
PID 113 0.813 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeInterceptAttType, 
      mdlStudent$lmeSlopesAttType) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 16: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudent\(lmeInterceptAttType </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 5 </td> <td style="text-align:center;"> 36684 </td> <td style="text-align:center;"> 36717 </td> <td style="text-align:center;"> -18337 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlStudent\)lmeSlopesAttType 2 10 36476 36542 -18228 1 vs 2 217.9 0 
# Save variances
mdlStudent$varSlopesAttType <- 
  lme4::VarCorr(mdlStudent$lmeSlopesAttType)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwSlopesAttType <- 
  1 - (as.numeric(mdlStudent$varSlopesAttType[1]) / as.numeric(mdlStudent$varInterceptAttType[1]))
mdlStudent$varBtwPercSlopesAttType <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttType[1]) / as.numeric(mdlStudent$varInterceptAttType[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinSlopesAttType <- 
  1 - (as.numeric(mdlStudent$varSlopesAttType[2]) / as.numeric(mdlStudent$varInterceptAttType[2]))
mdlStudent$varWithinPercSlopesAttType <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttType[2]) / as.numeric(mdlStudent$varInterceptAttType[2]))) * 100

# Assumption Checks:
mdlStudent$diagSlopesAttType <-
  sjPlot::plot_model(mdlStudent$lmerSlopesAttType, type = "diag")
grid.arrange(
  mdlStudent$diagSlopesAttType[[1]],
  mdlStudent$diagSlopesAttType[[2]]$`PID`,
  mdlStudent$diagSlopesAttType[[3]],
  mdlStudent$diagSlopesAttType[[4]]
)

# Plot prediction model
mdlStudent$predSlopesAttType <- 
  studentInteractionType %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudent$lmeSlopesAttType,
                           studentInteractionType,
                           re.form = NA
                           )
         )

(
  mdlStudent$predPltSlopesAttType <-
    ggplot(data = mdlStudent$predSlopesAttType %>% filter(PID %in% studentPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Student_PredictionPlot_SlopesAttType.png",
  mdlStudent$predPltSlopesAttType,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

TL;DR: Random slopes don’t add much for this super simple model.

Interaction Frequency and Interaction Quality

\[\begin{equation} \tag{6} Attitude = ContactFreq \times AverageContactQual \end{equation}\]

# correlation panel
pairs.panels.new(
  studentContactFreq %>% select(SumContactNL, SumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)",
    "Sum:\nNumber of Outgroup Contacts (NL)",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

# correlation panel with interaction sums winsorized
pairs.panels.new(
  studentContactFreq %>% select(WinSumContactNL, WinSumContactNLAll, AvQuality, AvAttitude),
  labels = c(
    "Sum:\nNumer of beeps with Outgroup Contact (NL)\n[Winsorized]",
    "Sum:\nNumber of Outgroup Contacts (NL)\n[Winsorized]",
    "Mean:\nInteraction Quality",
    "Mean:\nOutgroup Attitudes (NL)"
  )
)

TL;DR: There is no meaningfull correlation between the participants’ Average Interaction Quality and their Average Outgroup Attitudes. However, the frequency of intergroup interactions had a significant correlation with both the average interaction quality and average outgroup attitudes.

# regression
mdlStudent$lmAttFreqQualX <-
  lm(AvAttitude ~ SumContactNL * AvQuality, data = studentContactFreq)
# summary(mdlStudent$lmAttFreqQualX)

summ(
  mdlStudent$lmAttFreqQualX,
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 113
Dependent variable AvAttitude
Type OLS linear regression
F(3,109) 3.748
R² 0.094
Adj. R² 0.069
Est. 2.5% 97.5% t val. p
(Intercept) 67.618 64.156 71.080 38.707 0.000
SumContactNL 0.805 0.320 1.290 3.291 0.001
AvQuality 0.253 -0.135 0.641 1.292 0.199
SumContactNL:AvQuality 0.023 -0.032 0.078 0.819 0.414
Standard errors: OLS; Continuous predictors are mean-centered. 
interactions::interact_plot(
  mdlStudent$lmAttFreqQualX,
  pred = SumContactNL,
  modx = AvQuality,
  interval = TRUE,
  partial.residuals = TRUE
)

interactions::johnson_neyman(mdlStudent$lmAttFreqQualX,
                             pred = SumContactNL,
                             modx = AvQuality,
                             alpha = .05)
## JOHNSON-NEYMAN INTERVAL 
## 
## When AvQuality is INSIDE the interval [73.80, 100.18], the slope of SumContactNL is p < .05.
## 
## Note: The range of observed values of AvQuality is [61.50, 100.00]

TL;DR: We find that, although the effect goes into the right direction (i.e., higher average interaction quality increases the positive effect of intergroup contact), there was no effect of average interaction quality (neither directly nor in a moderating role of interaction frequency). However, this regression ignores within person relationships of these variables (i.e., summarizes all within person variances into average scores).

Need Fulfillment

Need fulfillment and Attitudes

Unconditional model

make a model with only the data but a random intercept (unconditional model) for the outgroup interaction subset.

# see how large our outgroup interaction subset actually is
tbl_cross(
  studentInteractionType,
  row = OutgroupInteraction,
  col = NonOutgroupInteraction,
  percent = "row"
)
Characteristic NonOutgroupInteraction Total
No Yes
OutgroupInteraction
No 1,695 (42%) 2,335 (58%) 4,030 (100%)
Yes 329 (35%) 606 (65%) 935 (100%)
Total 2,024 (41%) 2,941 (59%) 4,965 (100%)
 
# create outgroup interaction subset
studentOutgroupInteraction <- studentInteractionType %>%
  filter(OutgroupInteraction == "Yes")

# create empty list to organize models
mdlStudentOut <- 
  list(
    Att = list(),
    Qlt = list()
  )

\[\begin{equation} \tag{23} Attitude_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmerNull <-
  lme4::lmer(AttitudesDutch ~ 1 + (1 | PID), 
             data = studentOutgroupInteraction) # use optim if it does not converge
mdlStudentOut$Att$lmeNull <-
  lme(
    AttitudesDutch ~ 1,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.ML.r) #or with the lme function
summ(mdlStudentOut$Att$lmerNull, digits = 3, center = TRUE)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7250.070
BIC 7264.591
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.724
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 70.736 1.609 43.954 102.806 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 16.042
Residual 9.898
Grouping Variables
Group # groups ICC
PID 108 0.724 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudentOut$Att$varNull <- 
  VarCorr(mdlStudentOut$Att$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudentOut$Att$tauNull <- 
  as.numeric(mdlStudentOut$Att$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudentOut$Att$sigmaNull <- 
  as.numeric(mdlStudentOut$Att$varNull[2])
# The ICC estimate (between/between+within) is:
mdlStudentOut$Att$IccNull <-
  (as.numeric(mdlStudentOut$Att$varNull[1]) / (as.numeric(mdlStudentOut$Att$varNull[1]) + as.numeric(mdlStudentOut$Att$varNull[2])))
mdlStudentOut$Att$IccPercNull <-
  ((as.numeric(mdlStudentOut$Att$varNull[1]) / (as.numeric(mdlStudentOut$Att$varNull[1]) + as.numeric(mdlStudentOut$Att$varNull[2])))) * 100

random intercept

\[\begin{equation} \tag{8} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCore <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(AttitudesDutch ~ KeyNeedFullfillment + (1 | PID), 
       data = studentOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7234.687
BIC 7254.049
Pseudo-R² (fixed effects) 0.009
Pseudo-R² (total) 0.728
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.811 67.685 73.936 44.404 102.944 0.000
KeyNeedFullfillment 0.096 0.057 0.134 4.866 851.477 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.900
Residual 9.777
Grouping Variables
Group # groups ICC
PID 108 0.726 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeNull, 
      mdlStudentOut$Att$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 17: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622 NA NA
mdlStudentOut\(Att\)lmeInterceptCore 2 4 7235 7254 -3613 1 vs 2 17.38 0 
# Save variances
mdlStudentOut$Att$varInterceptCore <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCore <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))
mdlStudentOut$Att$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCore <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))
mdlStudentOut$Att$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCore[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))) * 100

TL;DR: Situational core need relates significantly to outgroup attitudes. The core need has little explained variance between participants (compared to the null model: Variance Explained = 1 – (Var with Predictor/Var without Predictor); 1.77%). The variance explained within participants is small to medium (2.44%)

random slope

\[\begin{equation} \tag{9} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Att$lmeSlopesCore <-
  lme(
    AttitudesDutch ~
      KeyNeedFullfillment,
    random = ~ 1 + KeyNeedFullfillment | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCore <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment +
      (1 + KeyNeedFullfillment | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7218.921
BIC 7247.965
Pseudo-R² (fixed effects) 0.017
Pseudo-R² (total) 0.745
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.915 67.809 74.022 44.739 102.943 0.000
KeyNeedFullfillment 0.132 0.073 0.191 4.382 45.436 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.774
PID KeyNeedFullfillment 0.150
Residual 9.472
Grouping Variables
Group # groups ICC
PID 108 0.735 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeInterceptCore, 
      mdlStudentOut$Att$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 18: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeInterceptCore 1 4 7235 7254 -3613 NA NA
mdlStudentOut\(Att\)lmeSlopesCore 2 6 7219 7248 -3603 1 vs 2 19.76 0 
# Save variances
mdlStudentOut$Att$varSlopesCore <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))
mdlStudentOut$Att$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))
mdlStudentOut$Att$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCore[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCore <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCore[[1]],
  mdlStudentOut$Att$diagSlopesCore[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCore[[3]],
  mdlStudentOut$Att$diagSlopesCore[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCore <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCore,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCore <-
    ggplot(data = mdlStudentOut$Att$predSlopesCore %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCore.png",
  mdlStudentOut$Att$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

TL;DR: The random slope adds significantly to the prediction model.

Need fulfillment and Interaction Quality

Unconditional model

make a model with only the data but a random intercept (unconditional model) for the outgroup interaction subset.

\[\begin{equation} \tag{24} InteractionQuality_{ti} = \gamma_{00} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Qlt$lmerNull <-
  lme4::lmer(quality_overall ~ 1 + (1 | PID), 
             data = studentOutgroupInteraction) # use optim if it does not converge
mdlStudentOut$Qlt$lmeNull <-
  lme(
    quality_overall ~ 1,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    control = list(opt = "nlmimb")
  ) # use optim if it does not converge

# Get summary with p-values (Satterthwaite's method)
# summary(Null.Out.Qual.ML.r) #or with the lme function
summ(mdlStudentOut$Qlt$lmerNull, digits = 3, center = TRUE)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8179.693
BIC 8194.215
Pseudo-R² (fixed effects) 0.000
Pseudo-R² (total) 0.144
Fixed Effects
Est. S.E. t val. d.f. p
(Intercept) 78.849 1.017 77.536 100.065 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 7.506
Residual 18.318
Grouping Variables
Group # groups ICC
PID 108 0.144 
# generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch~1 + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000,
#                  parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-Null-CI.csv")

# Save variances
mdlStudentOut$Qlt$varNull <- 
  VarCorr(mdlStudentOut$Qlt$lmeNull) # save variances
# The estimate of (between-group or Intercept variance, tau_{00}^2):
mdlStudentOut$Qlt$tauNull <- 
  as.numeric(mdlStudentOut$Qlt$varNull[1])
# and the estimate of (within-group or residual variance, sigma^2) is:
mdlStudentOut$Qlt$sigmaNull <- 
  as.numeric(mdlStudentOut$Qlt$varNull[2])
# The ICC estimate (between/between+within) is:
mdlStudentOut$Qlt$IccNull <-
  (as.numeric(mdlStudentOut$Qlt$varNull[1]) / (as.numeric(mdlStudentOut$Qlt$varNull[1]) + as.numeric(mdlStudentOut$Qlt$varNull[2])))
mdlStudentOut$Qlt$IccPercNull <-
  ((as.numeric(mdlStudentOut$Qlt$varNull[1]) / (as.numeric(mdlStudentOut$Qlt$varNull[1]) + as.numeric(mdlStudentOut$Qlt$varNull[2])))) * 100

random intercept

\[\begin{equation} \tag{11} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Qlt$lmeInterceptCore <-
  lme(
    quality_overall ~ KeyNeedFullfillment,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudentOut$Qlt$lmerInterceptCore <- lmer(quality_overall ~ KeyNeedFullfillment + (1 | PID), 
       data = studentOutgroupInteraction),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8081.658
BIC 8101.021
Pseudo-R² (fixed effects) 0.106
Pseudo-R² (total) 0.221
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 78.915 77.091 80.739 84.803 102.614 0.000
KeyNeedFullfillment 0.348 0.284 0.413 10.550 932.928 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.678
Residual 17.389
Grouping Variables
Group # groups ICC
PID 108 0.129 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(quality_overall ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Qlt$lmeNull, 
      mdlStudentOut$Qlt$lmeInterceptCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 19: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Qlt\)lmeNull 1 3 8180 8194 -4087 NA NA
mdlStudentOut\(Qlt\)lmeInterceptCore 2 4 8082 8101 -4037 1 vs 2 100 0 
# Save variances
mdlStudentOut$Qlt$varInterceptCore <-
  lme4::VarCorr(mdlStudentOut$Qlt$lmeInterceptCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Qlt$varBtwInterceptCore <- 
  1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]) / as.numeric(mdlStudentOut$Qlt$varNull[1]))
mdlStudentOut$Qlt$varBtwPercInterceptCore <- 
  (1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]) / as.numeric(mdlStudentOut$Qlt$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Qlt$varWithinInterceptCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]) / as.numeric(mdlStudentOut$Qlt$varNull[2]))
mdlStudentOut$Qlt$varWithinPercInterceptCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]) / as.numeric(mdlStudentOut$Qlt$varNull[2]))) * 100

TL;DR: Situational core need relates significantly to interaction quality. The core need explained substantial variance between participants (20.84%). The variance explained within participants is also medium (9.88%)

random slope

\[\begin{equation} \tag{12} InteractionQuality_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Qlt$lmeSlopesCore <-
  lme(
    quality_overall ~
      KeyNeedFullfillment,
    random = ~ 1 + KeyNeedFullfillment | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Qlt$lmerSlopesCore <-
    lmer(
      quality_overall ~
        KeyNeedFullfillment +
        (1 + KeyNeedFullfillment | PID),
      data = studentOutgroupInteraction
    ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable quality_overall
Type Mixed effects linear regression
AIC 8047.231
BIC 8076.275
Pseudo-R² (fixed effects) 0.143
Pseudo-R² (total) 0.335
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 78.882 77.114 80.651 87.409 102.262 0.000
KeyNeedFullfillment 0.416 0.306 0.526 7.400 52.447 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 6.080
PID KeyNeedFullfillment 0.351
Residual 16.525
Grouping Variables
Group # groups ICC
PID 108 0.119 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(mdlStudentOut$Qlt$lmerSlopesCore,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Qlt$lmeInterceptCore, 
      mdlStudentOut$Qlt$lmeSlopesCore) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 20: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Qlt\)lmeInterceptCore 1 4 8082 8101 -4037 NA NA
mdlStudentOut\(Qlt\)lmeSlopesCore 2 6 8047 8076 -4018 1 vs 2 38.42 0 
# Save variances
mdlStudentOut$Qlt$varSlopesCore <- 
  lme4::VarCorr(mdlStudentOut$Qlt$lmeSlopesCore)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Qlt$varBtwSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[1]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]))
mdlStudentOut$Qlt$varBtwPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[1]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Qlt$varWithinSlopesCore <-
  1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[2]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]))
mdlStudentOut$Qlt$varWithinPercSlopesCore <-
  (1 - (as.numeric(mdlStudentOut$Qlt$varSlopesCore[2]) / as.numeric(mdlStudentOut$Qlt$varInterceptCore[2]))) * 100

# Assumption Checks:
mdlStudentOut$Qlt$diagSlopesCore <-
  sjPlot::plot_model(mdlStudentOut$Qlt$lmerSlopesCore, type = "diag")
grid.arrange(
  mdlStudentOut$Qlt$diagSlopesCore[[1]],
  mdlStudentOut$Qlt$diagSlopesCore[[2]]$`PID`,
  mdlStudentOut$Qlt$diagSlopesCore[[3]],
  mdlStudentOut$Qlt$diagSlopesCore[[4]]
)

# Plot prediction model
mdlStudentOut$Qlt$predSlopesCore <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Qlt$lmeSlopesCore,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Qlt$predPltSlopesCore <-
    ggplot(data = mdlStudentOut$Qlt$predSlopesCore %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesCore.png",
  mdlStudentOut$Qlt$predPltSlopesCore,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

TL;DR: random slopes does add significantly.

Need fulfillment, Interaction Quality, and Attitudes

random intercept

\[\begin{equation} \tag{13} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCoreQlt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment + quality_overall,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  lmer(
    AttitudesDutch ~ KeyNeedFullfillment + quality_overall + (1 | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7179.308
BIC 7203.511
Pseudo-R² (fixed effects) 0.032
Pseudo-R² (total) 0.747
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.720 67.611 73.829 44.581 102.807 0.000
KeyNeedFullfillment 0.045 0.006 0.084 2.251 846.076 0.025
quality_overall 0.151 0.115 0.187 8.114 841.537 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.855
Residual 9.421
Grouping Variables
Group # groups ICC
PID 108 0.739 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeNull, 
  mdlStudentOut$Att$lmeInterceptCoreQlt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 21: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622 NA NA
mdlStudentOut\(Att\)lmeInterceptCoreQlt 2 5 7179 7203 -3585 1 vs 2 74.76 0 
# Save variances
mdlStudentOut$Att$varInterceptCoreQlt <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCoreQlt <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))
mdlStudentOut$Att$varBtwPercInterceptCoreQlt <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))
mdlStudentOut$Att$varWithinPercInterceptCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varNull[2]))) * 100

TL;DR: Variance explained in outgroup attitudes is of medium effect size (between: 2.32%, within: 9.42%)

random slope

\[\begin{equation} \tag{14} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}InteractionQuality_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudentOut$Att$lmeSlopesCoreQlt <-
  lme(
    AttitudesDutch ~
      KeyNeedFullfillment + quality_overall,
    random = ~ 1 + KeyNeedFullfillment + quality_overall | PID,
    control = lmeControl(opt = "optim"),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCoreQlt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment + quality_overall +
      (1 + KeyNeedFullfillment + quality_overall | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7137.672
BIC 7186.077
Pseudo-R² (fixed effects) 0.039
Pseudo-R² (total) 0.787
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.582 67.471 73.692 44.472 102.581 0.000
KeyNeedFullfillment 0.039 -0.004 0.082 1.789 24.796 0.086
quality_overall 0.173 0.114 0.233 5.723 50.585 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.844
PID KeyNeedFullfillment 0.055
PID quality_overall 0.192
Residual 8.733
Grouping Variables
Group # groups ICC
PID 108 0.767 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeNull,
  mdlStudentOut$Att$lmeInterceptCoreQlt,
  mdlStudentOut$Att$lmeSlopesCoreQlt
) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 22: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeNull 1 3 7250 7265 -3622 NA NA
mdlStudentOut\(Att\)lmeInterceptCoreQlt 2 5 7179 7203 -3585 1 vs 2 74.76 0
mdlStudentOut\(Att\)lmeSlopesCoreQlt 3 10 7138 7187 -3559 2 vs 3 51.11 0 
# Save variances
mdlStudentOut$Att$varSlopesCoreQlt <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCoreQlt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]))
mdlStudentOut$Att$varBtwPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCoreQlt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]))
mdlStudentOut$Att$varWithinPercSlopesCoreQlt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreQlt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreQlt[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCoreQlt <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCoreQlt, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCoreQlt[[1]],
  mdlStudentOut$Att$diagSlopesCoreQlt[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCoreQlt[[3]],
  mdlStudentOut$Att$diagSlopesCoreQlt[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCoreQlt <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCoreQlt,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCoreQlt <-
    ggplot(data = mdlStudentOut$Att$predSlopesCoreQlt %>% filter(PID %in% studentOutPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCoreQlt.png",
  mdlStudentOut$Att$predPltSlopesCoreQlt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

Compare models

Interpretation of Mediation like results.

Check for robustness

Interaction Type

random intercept

\[\begin{equation} \tag{15} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudent$lmeInterceptAttCoreInt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction,
    random =  ~ 1 | PID,
    data = studentInteractionType
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudent$lmerInterceptAttCoreInt <- lmer(
    AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction + (1 | PID),
    data = studentInteractionType
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36668.990
BIC 36708.051
Pseudo-R² (fixed effects) 0.004
Pseudo-R² (total) 0.802
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.772 63.373 70.171 38.501 112.066 0.000
KeyNeedFullfillment 0.025 0.012 0.039 3.587 4859.664 0.000
OutgroupInteraction 2.638 1.915 3.360 7.154 4861.666 0.000
KeyNeedFullfillment:OutgroupInteraction 0.049 0.014 0.085 2.742 4853.890 0.006
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.349
Residual 9.152
Grouping Variables
Group # groups ICC
PID 113 0.801 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentInteractionType),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 23: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudent\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 36741 </td> <td style="text-align:center;"> 36760 </td> <td style="text-align:center;"> -18367 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlStudent\)lmeInterceptAttCoreInt 2 6 36669 36708 -18328 1 vs 2 77.94 0 
# Save variances
mdlStudent$varInterceptAttCoreInt <- 
  lme4::VarCorr(mdlStudent$lmeInterceptAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwInterceptAttCoreInt <-
  1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[1]) / as.numeric(mdlStudent$varAttNull[1]))
mdlStudent$varBtwPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[1]) / as.numeric(mdlStudent$varAttNull[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinInterceptAttCoreInt <-
  1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[2]) / as.numeric(mdlStudent$varAttNull[2]))
mdlStudent$varWithinPercInterceptAttCoreInt <-
  (1 - (as.numeric(mdlStudent$varInterceptAttCoreInt[2]) / as.numeric(mdlStudent$varAttNull[2]))) * 100
random slope

\[\begin{equation} \tag{16} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}OutgroupInteraction_{ti} + \gamma_{30}KeyNeedFulfillXOutgroupInteraction_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model (optimizer needed to reach convergence)
mdlStudent$lmeSlopesAttCoreInt <- lme(
  AttitudesDutch ~
    KeyNeedFullfillment * OutgroupInteraction,
  random = ~ 1 + KeyNeedFullfillment + OutgroupInteraction | PID,
  control = lmeControl(opt = "optim"),
  data = studentInteractionType
)

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudent$lmerSlopesAttCoreInt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment * OutgroupInteraction +
      (1 + KeyNeedFullfillment + OutgroupInteraction | PID),
    data = studentInteractionType
  ), 
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 4965
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 36473.506
BIC 36545.117
Pseudo-R² (fixed effects) 0.006
Pseudo-R² (total) 0.815
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 66.773 63.299 70.248 37.666 113.093 0.000
KeyNeedFullfillment 0.029 0.011 0.047 3.166 84.092 0.002
OutgroupInteraction 2.898 1.365 4.431 3.705 94.703 0.000
KeyNeedFullfillment:OutgroupInteraction 0.060 0.021 0.099 3.039 1342.328 0.002
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 18.752
PID KeyNeedFullfillment 0.055
PID OutgroupInteraction 6.932
Residual 8.794
Grouping Variables
Group # groups ICC
PID 113 0.820 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudent$lmeAttNull, 
      mdlStudent$lmeInterceptAttCoreInt,
      mdlStudent$lmeSlopesAttCoreInt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 24: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudent\(lmeAttNull </td> <td style="text-align:center;"> 1 </td> <td style="text-align:center;"> 3 </td> <td style="text-align:center;"> 36741 </td> <td style="text-align:center;"> 36760 </td> <td style="text-align:center;"> -18367 </td> <td style="text-align:center;"> </td> <td style="text-align:center;"> NA </td> <td style="text-align:center;"> NA </td> </tr> <tr> <td style="text-align:left;"> mdlStudent\)lmeInterceptAttCoreInt 2 6 36669 36708 -18328 1 vs 2 77.94 0
mdlStudent$lmeSlopesAttCoreInt 3 11 36474 36545 -18226 2 vs 3 205.46 0 
# Save variances
mdlStudent$varSlopesAttCoreInt <- 
  lme4::VarCorr(mdlStudent$lmeSlopesAttCoreInt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudent$varBtwSlopesAttCoreInt <- 
  1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[1]) / as.numeric(mdlStudent$varInterceptAttCoreInt[1]))
mdlStudent$varBtwPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[1]) / as.numeric(mdlStudent$varInterceptAttCoreInt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudent$varWithinSlopesAttCoreInt <- 
  1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[2]) / as.numeric(mdlStudent$varInterceptAttCoreInt[2]))
mdlStudent$varWithinPercSlopesAttCoreInt <- 
  (1 - (as.numeric(mdlStudent$varSlopesAttCoreInt[2]) / as.numeric(mdlStudent$varInterceptAttCoreInt[2]))) * 100

# Assumption Checks:
mdlStudent$diagSlopesAttCoreInt <-
  sjPlot::plot_model(mdlStudent$lmerSlopesAttCoreInt, type = "diag")
grid.arrange(
  mdlStudent$diagSlopesAttCoreInt[[1]],
  mdlStudent$diagSlopesAttCoreInt[[2]]$`PID`,
  mdlStudent$diagSlopesAttCoreInt[[3]],
  mdlStudent$diagSlopesAttCoreInt[[4]]
)

# Plot prediction model
mdlStudent$predSlopesAttCoreInt <- 
  studentInteractionType %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudent$lmeSlopesAttCoreInt,
                           studentInteractionType,
                           re.form = NA
                           )
         )

(
  mdlStudent$predPltSlopesAttCoreInt <-
    ggplot(data = mdlStudent$predSlopesAttCoreInt %>% filter(PID %in% studentPltIDs), 
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap(~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/Student_PredictionPlot_SlopesAttCoreInt.png",
  mdlStudent$predPltSlopesAttCoreInt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)
Plot Interaction
# visualize interaction
## Without ML structure
ggplot(data = studentInteractionType,
       aes(x = KeyNeedFullfillment,
           y = AttitudesDutch,
           fill = OutgroupInteraction)) +
  #geom_point()+
  geom_smooth(method = 'lm',
              aes(linetype = OutgroupInteraction),
              color = "black") +
  #facet_wrap(~PID, ncol = 6)+
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  #scale_colour_manual(values=c("grey20", "black"), name="Intergroup Contact")+
  scale_y_continuous(
    limits = c(50, 100),
    breaks = seq(50, 100, by = 10),
    position = "left"
  ) +
  scale_x_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(
    title = "Without ML stucture",
    x = "Fulfillment Core Need",
    y = "Outgroup Attitudes",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact"
  ) +
  theme_Publication() +
  theme(legend.position = "bottom",
        legend.key.size = unit(1, "cm"))

## With ML structure
# create parameters for prediction
datNew = data.frame(
  OutgroupInteraction = factor(rep(c("No", "Yes"), each = 21)),
  KeyNeedFullfillment = rep(seq(0, 100, 5), 2),
  PID = 0
)

# Predict values, clean up and calculate SE
PI <-
  merTools::predictInterval(
    merMod = mdlStudent$lmerSlopesAttCoreInt,
    newdata = datNew,
    level = 0.95,
    stat = "mean",
    type = "linear.prediction",
    include.resid.var = F,
    fix.intercept.variance = F
  )
mdlStudent$predInterceptAttCoreIntX <- 
  cbind(datNew, PI)
mdlStudent$predInterceptAttCoreIntX$se <-
  (mdlStudent$predInterceptAttCoreIntX$upr - mdlStudent$predInterceptAttCoreIntX$fit) / 1.96
rm(datNew, PI)

# Plot predicted values with SE
ggplot(
  mdlStudent$predInterceptAttCoreIntX,
  aes(x = KeyNeedFullfillment,
      y = fit,
      fill = OutgroupInteraction)
)+
  #geom_point() +
  geom_line(aes(linetype = as.factor(OutgroupInteraction)), size = 1) +
  #facet_wrap(~PID, ncol = 6)+
  geom_ribbon(data = mdlStudent$predInterceptAttCoreIntX,
              aes(ymin = fit - se, ymax = fit + se),
              alpha = 0.3) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
  scale_linetype_manual(values = c("dashed", "solid")) +
  scale_fill_manual(values = c("darkgrey", "black")) +
  labs(
    x = "Fulfillment Core Need",
    y = "Outgroup Attitude (NL)",
    fill = "Intergroup Contact",
    linetype = "Intergroup Contact",
    title = "Based on Model Predictions"
  ) +
  theme_Publication()

# #### Bayesian estimation !! ONLY RUN ON FINAL RENDER !! Takes forever ####
# options(mc.cores = parallel::detectCores())  # Run many chains simultaneously
# brmfit <- brm(
#   AttitudesDutch ~ KeyNeedFullfillment * OutgroupInteraction + 
#     (1 + KeyNeedFullfillment + OutgroupInteraction | PID),
#   data = studentInteractionType,
#   family = gaussian,
#   iter = 1000,
#   chains = 4
# )
# 
# # create parameters for prediction:
# datNew = data.frame(
#   OutgroupInteraction = rep(c("No", "Yes"), each = 101),
#   KeyNeedFullfillment = rep(seq(0, 100, 1), 2)
# )
# # Save predicted values and adjust names and labels
# fitavg <-
#   cbind(datNew,
#         fitted(brmfit, newdata = datNew, re_formula = NA)[, -2])
# names(fitavg)[names(fitavg) == "Estimate"] = "pred"
# fitavg$se <- (fitavg$Q97.5 - fitavg$pred) / 1.96
# 
# # Plot Bayesian SE prediction interval
# ggplot(fitavg,
#        aes(x = KeyNeedFullfillment,
#            y = pred,
#            fill = OutgroupInteraction)) +
#   scale_x_continuous(breaks = seq(0, 100, 10)) +
#   scale_y_continuous(limits = c(50, 100), breaks = seq(50, 100, 10)) +
#   geom_line(aes(linetype = OutgroupInteraction), size = 1) +
#   geom_ribbon(aes(ymin = pred - se, ymax = pred + se), alpha = 0.2) +
#   scale_linetype_manual(values = c("dashed", "solid")) +
#   scale_fill_manual(values = c("darkgrey", "black")) +
#   labs(
#     x = "Fulfillment Core Need",
#     y = "Outgroup Attitude (NL)",
#     fill = "Intergroup Contact",
#     linetype = "Intergroup Contact",
#     title = "Based on Bayesian Prediction Interval"
#   ) +
#   theme_Publication()
# 
# # plot all overlayed posteriors:
# pst <- posterior_samples(brmfit, "b")
# ggplot(studentInteractionType,
#        aes(x = KeyNeedFullfillment, y = AttitudesDutch)) +
#   geom_point(shape = 4, alpha = .1) +
#   geom_tile() +
#   geom_abline(
#     data = pst,
#     aes(intercept = b_Intercept, slope = b_KeyNeedFullfillment),
#     alpha = .025,
#     size = .4
#   ) +
#   labs(title = "slope Posteriors",
#        x = "Fulfillment Core Need",
#        y = "Outgroup Attitudes (NL)") +
#   theme_Publication()
# rm(datNew, brmfit, fitavg, pst)

Other psychological needs

random intercept

\[\begin{equation} \tag{17} Attitude_{ti} = \gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + e_{ti} \\ \textrm{with}\ u_{0i} \sim \mathcal{N}(0,\tau_{00}^2)\ \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Create and save Model
mdlStudentOut$Att$lmeInterceptCoreSdt <-
  lme(
    AttitudesDutch ~ KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction,
    random = ~ 1 | PID,
    data = studentOutgroupInteraction,
    na.action = na.exclude
  )

# Get summary with p-values (Satterthwaite's method)
summ(
  mdlStudentOut$Att$lmerInterceptCoreSdt <- lmer(
    AttitudesDutch ~ KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction + (1 | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7212.060
BIC 7245.944
Pseudo-R² (fixed effects) 0.028
Pseudo-R² (total) 0.745
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.753 67.624 73.881 44.325 102.929 0.000
KeyNeedFullfillment 0.062 0.023 0.101 3.147 845.995 0.002
Competence 0.056 0.013 0.098 2.575 852.148 0.010
Autonomy 0.034 -0.015 0.083 1.344 870.215 0.179
RelatednessInteraction 0.058 0.033 0.084 4.425 848.412 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.943
Residual 9.519
Grouping Variables
Group # groups ICC
PID 108 0.737 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(lmer(AttitudesDutch ~ Contact_dum + inNonDutch + (1|PID),data=studentOutgroupInteraction),
#                  method="boot",nsim=1000, parallel = "multicore",
#                  ncpus = 4, seed = 42),
#          "output/tables/ML-Inter-CI.csv")

# Compare new model to previous step
anova(
  mdlStudentOut$Att$lmeInterceptCore, 
  mdlStudentOut$Att$lmeInterceptCoreSdt
  ) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 25: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeInterceptCore 1 4 7235 7254 -3613 NA NA
mdlStudentOut\(Att\)lmeInterceptCoreSdt 2 7 7212 7246 -3599 1 vs 2 28.63 0 
# Save variances
mdlStudentOut$Att$varInterceptCoreSdt <-
  lme4::VarCorr(mdlStudentOut$Att$lmeInterceptCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwInterceptCoreSdt <- 
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))
mdlStudentOut$Att$varBtwPercInterceptCoreSdt <- 
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCore[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinInterceptCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))
mdlStudentOut$Att$varWithinPercInterceptCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCore[2]))) * 100
random slope

\[\begin{equation} \tag{18} Attitude_{ti} = \\gamma_{00} + \gamma_{10}KeyNeedFulfill_{ti} + \gamma_{20}Autonomy_{ti} + \gamma_{30}Competence_{ti} + \gamma_{40}Relatedness_{ti} + u_{0i} + u_{1i} + e_{ti} \\ \textrm{with}\ \begin{bmatrix} u_{0i}\\ u_{1i}\end{bmatrix} \sim \mathcal{N}(0,\Omega_u): \Omega_u=\begin{bmatrix} \tau_{0}^2 & \\ \tau_{01} & \tau_{1}^2 \end{bmatrix} \textrm{and}\ e_{ti} \sim \mathcal{N}(0,\sigma^2) \end{equation}\]

# Check Correlations 
studentOutgroupInteraction %>% 
  select(AttitudesDutch, KeyNeedFullfillment, Competence, Autonomy, RelatednessInteraction) %>%
  pairs.panels.new

studentRed <- 
  studentOutgroupInteraction %>% 
  # select(PID, AttitudesDutch, KeyNeedFullfillment, Competence, Autonomy, RelatednessInteraction) %>%
  group_by(PID) %>%
  filter(
    sd(AttitudesDutch) != 0,
    sd(KeyNeedFullfillment) != 0,
    sd(Competence) != 0,
    sd(Autonomy) != 0,
    sd(RelatednessInteraction) != 0
  ) %>%
  ungroup

studentRedCor <- 
  studentOutgroupInteraction %>%
  group_by(PID) %>%
  summarise(
    rAttCore = cor(AttitudesDutch, KeyNeedFullfillment),
    rAttComp = cor(AttitudesDutch, Competence),
    rAttAut = cor(AttitudesDutch, Autonomy),
    rAttRel = cor(AttitudesDutch, RelatednessInteraction),
    rCoreComp = cor(KeyNeedFullfillment, Competence),
    rCoreAut = cor(KeyNeedFullfillment, Autonomy),
    rCoreRel = cor(KeyNeedFullfillment, RelatednessInteraction),
    rCompAut = cor(Competence, Autonomy),
    rCompRel = cor(Competence, RelatednessInteraction),
    rAutRel = cor(Autonomy, RelatednessInteraction)
  ) %>%
  filter_at(vars(-PID), all_vars(abs(.) != "1"))
  # mutate(corMean = rowMeans(abs(.[2:ncol(.)]))) %>%
  # filter(corMean != "1")

studentRed2 <- 
  studentOutgroupInteraction %>%
  filter(PID %in% studentRedCor$PID)

# Create and save Model (optimizer needed to reach convergence) 
# Problem because some PPTs have 100 on all measures (SD = 0) AND/OR
# For some all cor = 1 or -1
# Removing these PPTs or their measurement beeps doesn't help
# However, removing eithe the Core Need, Autonomy, or Relatedness from the random slopes lets the model converge
# --> looks like there is an overfitting issue? But why does it work in lmer?
# FOR NOW: Autonomy removed from random slopes because weakest predictor
mdlStudentOut$Att$lmeSlopesCoreSdt <-
  nlme::lme(
    AttitudesDutch ~
      KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction,
    random = ~ 1 + KeyNeedFullfillment + Competence + RelatednessInteraction | PID, # Autonomy + 
    #method = "ML",
    #control = lmeControl(opt = "optim", optimMethod = "L-BFGS-B"),
    control = lmeControl(opt = "optim", maxIter = 100, msMaxIter = 100),
    data = studentOutgroupInteraction
  )

# Get summary with p-values (Satterthwaite's method) [+ save model for plotting]
summ(
  mdlStudentOut$Att$lmerSlopesCoreSdt <- lmer(
    AttitudesDutch ~
      KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction +
      (1 + KeyNeedFullfillment + Competence + Autonomy + RelatednessInteraction | PID),
    data = studentOutgroupInteraction
  ),
  confint = TRUE,
  digits = 3,
  center = TRUE
)
Observations 935
Dependent variable AttitudesDutch
Type Mixed effects linear regression
AIC 7214.338
BIC 7315.990
Pseudo-R² (fixed effects) 0.034
Pseudo-R² (total) 0.765
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 70.672 67.578 73.766 44.775 103.278 0.000
KeyNeedFullfillment 0.085 0.031 0.139 3.072 40.508 0.004
Competence 0.052 0.008 0.096 2.303 92.107 0.024
Autonomy 0.027 -0.022 0.076 1.083 254.850 0.280
RelatednessInteraction 0.066 0.035 0.097 4.180 40.814 0.000
p values calculated using Satterthwaite d.f. ; Continuous predictors are mean-centered.
Random Effects
Group Parameter Std. Dev.
PID (Intercept) 15.703
PID KeyNeedFullfillment 0.124
PID Competence 0.047
PID Autonomy 0.022
PID RelatednessInteraction 0.069
Residual 9.089
Grouping Variables
Group # groups ICC
PID 108 0.749 
# Generate 95% parametric bootstrap CIs (and save them as a csv-file):
# write.csv(confint(model.ran0,
#               method="boot",nsim=1000,
#               parallel = "multicore", ncpus = 4, seed = 42),
#          "output/tables/ML-RandomSlopes-CI.csv")

# Compare new model to previous step
anova(mdlStudentOut$Att$lmeInterceptCoreSdt, 
      mdlStudentOut$Att$lmeSlopesCoreSdt) %>%
  as.data.frame() %>%
  select(-call) %>%
  kbl(
    .,
    caption = "Student: Model Comparison",
    format = "html",
    linesep = "",
    booktabs = T,
    align = rep("c", ncol(.)),
    digits = 3
  ) %>%
  kable_styling(position = "left")
Table 26: Student: Model Comparison
Model df AIC BIC logLik Test L.Ratio p-value
mdlStudentOut\(Att\)lmeInterceptCoreSdt 1 7 7212 7246 -3599 NA NA
mdlStudentOut\(Att\)lmeSlopesCoreSdt 2 16 7207 7284 -3587 1 vs 2 23.1 0.006 
# Save variances
mdlStudentOut$Att$varSlopesCoreSdt <- 
  lme4::VarCorr(mdlStudentOut$Att$lmeSlopesCoreSdt)

# The estimate of between-group (or Intercept variance) explained:
# Variance Explained = 1 – (Var with Predictor/Var without Predictor)
mdlStudentOut$Att$varBtwSlopesCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]))
mdlStudentOut$Att$varBtwPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[1]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[1]))) * 100
# and the estimate of within-group (or residual variance) explained is:
mdlStudentOut$Att$varWithinSlopesCoreSdt <-
  1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]))
mdlStudentOut$Att$varWithinPercSlopesCoreSdt <-
  (1 - (as.numeric(mdlStudentOut$Att$varSlopesCoreSdt[2]) / as.numeric(mdlStudentOut$Att$varInterceptCoreSdt[2]))) * 100

# Assumption Checks:
mdlStudentOut$Att$diagSlopesCoreSdt <- 
  sjPlot::plot_model(mdlStudentOut$Att$lmerSlopesCoreSdt, type = "diag")
grid.arrange(
  mdlStudentOut$Att$diagSlopesCoreSdt[[1]],
  mdlStudentOut$Att$diagSlopesCoreSdt[[2]]$`PID`,
  mdlStudentOut$Att$diagSlopesCoreSdt[[3]],
  mdlStudentOut$Att$diagSlopesCoreSdt[[4]]
)

# Plot prediction model
mdlStudentOut$Att$predSlopesCoreSdt <- 
  studentOutgroupInteraction %>%
  select(AttitudesDutch, TIDnum, PID) %>% 
  mutate(measure = predict(mdlStudentOut$Att$lmeSlopesCoreSdt,
                           studentOutgroupInteraction,
                           re.form = NA
                           )
         )

(
  mdlStudentOut$Att$predPltSlopesCoreSdt <-
    ggplot(data = mdlStudentOut$Att$predSlopesCoreSdt %>% filter(PID %in% studentOutPltIDs),
           aes(x = TIDnum, y = measure)) +
    geom_line(alpha = 1, color = "blue") +
    geom_line(aes(y = AttitudesDutch), alpha = 1) +
    facet_wrap( ~ PID, ncol = 6) +
    xlab("Time") +
    ylab("Outgroup Attitudes") +
    theme_Publication()
)

ggsave(
  filename = "Figures/StudentOut_PredictionPlot_SlopesAttCoreStd.png",
  mdlStudentOut$Att$predPltSlopesCoreSdt,
  width = 18,
  height = 12,
  dpi = 800,
  units = "cm",
  device = "png"
)

Young Medical Professional Sample

Contact Hypothesis

Allport’s Conditions

Need Fulfillment

Software Information

The full session information with all relevant system information and all loaded and installed packages is available in the collapsible section below.

System Info
Table 27: R environment session info for reproducibility of results
Setting Value
version R version 4.1.1 (2021-08-10)
os macOS Big Sur 10.16
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Amsterdam
date 2021-10-18

Package Info
Table 28: Package info for reproducibility of results
Package Loaded version Date Source
bookdown 0.24 2021-09-02 CRAN (R 4.1.0)
brms 2.16.1 2021-08-23 CRAN (R 4.1.0)
data.table 1.14.0 2021-02-21 CRAN (R 4.1.0)
devtools 2.4.2 2021-06-07 CRAN (R 4.1.0)
dplyr 1.0.7 2021-06-18 CRAN (R 4.1.0)
ellipse 0.4.2 2020-05-27 CRAN (R 4.1.0)
Formula 1.2-4 2020-10-16 CRAN (R 4.1.0)
ggpattern 0.2.0 2021-10-11 Github ()
ggplot2 3.3.5 2021-06-25 CRAN (R 4.1.0)
ggthemes 4.2.4 2021-01-20 CRAN (R 4.1.0)
gridExtra 2.3 2017-09-09 CRAN (R 4.1.0)
gtsummary 1.4.2 2021-07-13 CRAN (R 4.1.0)
haven 2.4.3 2021-08-04 CRAN (R 4.1.0)
Hmisc 4.5-0 2021-02-28 CRAN (R 4.1.0)
jtools 2.1.4 2021-09-03 CRAN (R 4.1.0)
kableExtra 1.3.4 2021-02-20 CRAN (R 4.1.0)
knitr 1.36 2021-09-29 CRAN (R 4.1.0)
lattice 0.20-44 2021-05-02 CRAN (R 4.1.1)
lme4 1.1-27.1 2021-06-22 CRAN (R 4.1.0)
lubridate 1.7.10 2021-02-26 CRAN (R 4.1.0)
mada 0.5.10 2020-05-25 CRAN (R 4.1.0)
Matrix 1.3-4 2021-06-01 CRAN (R 4.1.1)
mvmeta 1.0.3 2019-12-10 CRAN (R 4.1.0)
mvtnorm 1.1-2 2021-06-07 CRAN (R 4.1.0)
nlme 3.1-152 2021-02-04 CRAN (R 4.1.1)
pander 0.6.4 2021-06-13 CRAN (R 4.1.0)
papaja 0.1.0.9997 2021-10-11 Github ()
plotly 4.9.4.9000 2021-08-28 Github ()
plyr 1.8.6 2020-03-03 CRAN (R 4.1.0)
psych 2.1.6 2021-06-18 CRAN (R 4.1.0)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 4.1.0)
Rcpp 1.0.7 2021-07-07 CRAN (R 4.1.0)
remedy 0.1.0 2018-12-03 CRAN (R 4.1.0)
reshape2 1.4.4 2020-04-09 CRAN (R 4.1.0)
rmarkdown 2.11 2021-09-14 CRAN (R 4.1.1)
sessioninfo 1.1.1 2018-11-05 CRAN (R 4.1.0)
stringi 1.7.5 2021-10-04 CRAN (R 4.1.1)
stringr 1.4.0 2019-02-10 CRAN (R 4.1.0)
survival 3.2-12 2021-08-13 CRAN (R 4.1.1)
tibble 3.1.5 2021-09-30 CRAN (R 4.1.0)
tidyr 1.1.4 2021-09-27 CRAN (R 4.1.0)
usethis 2.0.1 2021-02-10 CRAN (R 4.1.0)

Full Session Info (including loaded but unattached packages — for troubleshooting only)

R version 4.1.1 (2021-08-10)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages:

  • grid
  • stats
  • graphics
  • grDevices
  • datasets
  • utils
  • methods
  • base

other attached packages:

  • lubridate(v.1.7.10)
  • reshape2(v.1.4.4)
  • stringi(v.1.7.5)
  • stringr(v.1.4.0)
  • papaja(v.0.1.0.9997)
  • kableExtra(v.1.3.4)
  • Hmisc(v.4.5-0)
  • Formula(v.1.2-4)
  • survival(v.3.2-12)
  • lattice(v.0.20-44)
  • tidyr(v.1.1.4)
  • dplyr(v.1.0.7)
  • plyr(v.1.8.6)
  • data.table(v.1.14.0)
  • mada(v.0.5.10)
  • mvmeta(v.1.0.3)
  • ellipse(v.0.4.2)
  • mvtnorm(v.1.1-2)
  • devtools(v.2.4.2)
  • usethis(v.2.0.1)
  • pander(v.0.6.4)
  • tibble(v.3.1.5)
  • sessioninfo(v.1.1.1)
  • gtsummary(v.1.4.2)
  • jtools(v.2.1.4)
  • nlme(v.3.1-152)
  • lme4(v.1.1-27.1)
  • Matrix(v.1.3-4)
  • ggpattern(v.0.2.0)
  • gridExtra(v.2.3)
  • plotly(v.4.9.4.9000)
  • RColorBrewer(v.1.1-2)
  • haven(v.2.4.3)
  • ggthemes(v.4.2.4)
  • ggplot2(v.3.3.5)
  • psych(v.2.1.6)
  • brms(v.2.16.1)
  • Rcpp(v.1.0.7)
  • bookdown(v.0.24)
  • remedy(v.0.1.0)
  • knitr(v.1.36)
  • rmarkdown(v.2.11)

loaded via a namespace (and not attached):

  • estimability(v.1.3)
  • coda(v.0.19-4)
  • dygraphs(v.1.1.1.6)
  • rpart(v.4.1-15)
  • inline(v.0.3.19)
  • generics(v.0.1.0)
  • callr(v.3.7.0)
  • commonmark(v.1.7)
  • proxy(v.0.4-26)
  • chron(v.2.3-56)
  • tzdb(v.0.1.2)
  • webshot(v.0.5.2)
  • xml2(v.1.3.2)
  • httpuv(v.1.6.3)
  • StanHeaders(v.2.21.0-7)
  • assertthat(v.0.2.1)
  • xfun(v.0.26)
  • hms(v.1.1.1)
  • jquerylib(v.0.1.4)
  • bayesplot(v.1.8.1)
  • evaluate(v.0.14)
  • promises(v.1.2.0.1)
  • fansi(v.0.5.0)
  • igraph(v.1.2.6)
  • DBI(v.1.1.1)
  • tmvnsim(v.1.0-2)
  • htmlwidgets(v.1.5.4)
  • tensorA(v.0.36.2)
  • stats4(v.4.1.1)
  • purrr(v.0.3.4)
  • ellipsis(v.0.3.2)
  • crosstalk(v.1.1.1)
  • backports(v.1.2.1)
  • V8(v.3.4.2)
  • insight(v.0.14.3)
  • markdown(v.1.1)
  • RcppParallel(v.5.1.4)
  • vctrs(v.0.3.8)
  • remotes(v.2.4.0)
  • sjlabelled(v.1.1.8)
  • abind(v.1.4-5)
  • cachem(v.1.0.6)
  • withr(v.2.4.2)
  • checkmate(v.2.0.0)
  • emmeans(v.1.6.3)
  • xts(v.0.12.1)
  • prettyunits(v.1.1.1)
  • mnormt(v.2.0.2)
  • svglite(v.2.0.0)
  • cluster(v.2.1.2)
  • lazyeval(v.0.2.2)
  • crayon(v.1.4.1)
  • labeling(v.0.4.2)
  • pkgconfig(v.2.0.3)
  • pkgload(v.1.2.1)
  • blme(v.1.0-5)
  • nnet(v.7.3-16)
  • rlang(v.0.4.11)
  • lifecycle(v.1.0.1)
  • miniUI(v.0.1.1.1)
  • colourpicker(v.1.1.0)
  • modelr(v.0.1.8)
  • distributional(v.0.2.2)
  • rprojroot(v.2.0.2)
  • matrixStats(v.0.60.1)
  • datawizard(v.0.2.0)
  • loo(v.2.4.1)
  • boot(v.1.3-28)
  • zoo(v.1.8-9)
  • base64enc(v.0.1-3)
  • gamm4(v.0.2-6)
  • ggridges(v.0.5.3)
  • processx(v.3.5.2)
  • png(v.0.1-7)
  • viridisLite(v.0.4.0)
  • parameters(v.0.14.0)
  • rootSolve(v.1.8.2.2)
  • arm(v.1.11-2)
  • readr(v.2.0.2)
  • jpeg(v.0.1-9)
  • shinystan(v.2.5.0)
  • ggeffects(v.1.1.1)
  • scales(v.1.1.1)
  • memoise(v.2.0.0)
  • magrittr(v.2.0.1)
  • threejs(v.0.3.3)
  • compiler(v.4.1.1)
  • rstantools(v.2.1.1)
  • cli(v.3.0.1)
  • lmerTest(v.3.1-3)
  • interactions(v.1.1.5)
  • TMB(v.1.7.21)
  • ps(v.1.6.0)
  • Brobdingnag(v.1.2-6)
  • htmlTable(v.2.2.1)
  • MASS(v.7.3-54)
  • mgcv(v.1.8-36)
  • tidyselect(v.1.1.1)
  • forcats(v.0.5.1)
  • mixmeta(v.1.1.3)
  • projpred(v.2.0.2)
  • highr(v.0.9)
  • yaml(v.2.2.1)
  • latticeExtra(v.0.6-29)
  • bridgesampling(v.1.1-2)
  • sass(v.0.4.0)
  • tools(v.4.1.1)
  • lmom(v.2.8)
  • parallel(v.4.1.1)
  • rstudioapi(v.0.13)
  • foreach(v.1.5.1)
  • foreign(v.0.8-81)
  • gld(v.2.6.2)
  • posterior(v.1.1.0)
  • farver(v.2.1.0)
  • sjPlot(v.2.8.9)
  • digest(v.0.6.28)
  • shiny(v.1.6.0)
  • broom(v.0.7.9.9001)
  • performance(v.0.7.3)
  • later(v.1.3.0)
  • httr(v.1.4.2)
  • rsconnect(v.0.8.24)
  • effectsize(v.0.4.5)
  • sjstats(v.0.18.1)
  • colorspace(v.2.0-2)
  • rvest(v.1.0.1)
  • fs(v.1.5.0)
  • splines(v.4.1.1)
  • rematch2(v.2.1.2)
  • expm(v.0.999-6)
  • Exact(v.3.0)
  • renv(v.0.14.0)
  • shinythemes(v.1.2.0)
  • systemfonts(v.1.0.2)
  • xtable(v.1.8-4)
  • jsonlite(v.1.7.2)
  • nloptr(v.1.2.2.2)
  • rstan(v.2.21.2)
  • testthat(v.3.0.4)
  • gt(v.0.3.1)
  • R6(v.2.5.1)
  • broom.mixed(v.0.2.7)
  • pillar(v.1.6.3)
  • htmltools(v.0.5.2)
  • mime(v.0.12)
  • glue(v.1.4.2)
  • fastmap(v.1.1.0)
  • minqa(v.1.2.4)
  • DT(v.0.19)
  • class(v.7.3-19)
  • codetools(v.0.2-18)
  • pkgbuild(v.1.2.0)
  • utf8(v.1.2.2)
  • bslib(v.0.3.0)
  • numDeriv(v.2016.8-1.1)
  • pbkrtest(v.0.5.1)
  • curl(v.4.3.2)
  • DescTools(v.0.99.43)
  • gtools(v.3.9.2)
  • shinyjs(v.2.0.0)
  • glmmTMB(v.1.1.2.3)
  • merTools(v.0.5.2)
  • desc(v.1.3.0)
  • munsell(v.0.5.0)
  • e1071(v.1.7-9)
  • iterators(v.1.0.13)
  • broom.helpers(v.1.4.0)
  • sjmisc(v.2.8.7)
  • gtable(v.0.3.0)
  • bayestestR(v.0.10.5)


References